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INVESTIGATION  OF  THE  HIGHER-ORDER  ELEMENTS 
IN  THE  SAMS0N2  CODE 
Abstract 

by  Steven  Scott  Miller,  M.S. 

Washington  State  University 
August  1986 

Chair:  Harold  C.  Sorensen 

The  objective  of  this  research  effort  was  to  determine  if  the  current 
finite  element  formulation  in  the  SAMS0N2  code  for  the  eight-node  quadrila¬ 
teral  higher-order  isoparametric  continuum  element  produces  consistent  and 
reliable  results.  Some  effort  was  also  performed  toward  investigating  the 
finite  element  formulations  for  the  five-node  and  the  six-node  triangular  iso¬ 
parametric  continuum  elements. 

Four  tasks  were  executed  in  this  research  effort.  The  first  task  was  a 
study  of  some  of  the  previous  work  which  had  been  performed.  The  second  task 
involved  comparisons  between  program  results  and  analytical  solutions  for 
various  problems.  The  third  task  was  the  verification  and  correction  where 
necessary  of  the  finite  element  formulation  for  the  eight-node  quadrilateral 
continuum  element.  The  fourth  task  was  a  study  of  the  effects  on  the  problem 
results  which  were  caused  by  the  corrected  finite  element  formulation. 

The  results  from  this  investigation  showed  that  the  current  finite  ele¬ 
ment  formulation  for  the  eight-node  quadrilateral  continuum  element  produces 
accurate  results  for  problems  which  are  analyzed  using  plane  continuum  ele¬ 
ments.  However,  corrections  to  the  finite  element  formulation  for  the  analy¬ 
ses  of  axi symmetric  conti nua  were  necessary  in  order  for  accurate  results  to 
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Introduction 


The  content  of  this  report  gives  a  detailed  description  of  the 
research  which  was  performed  during  the  investigation  of  the 
higher-order  isoparametric  continuum  finite  elements  (five-node 
triangle,  six-node  triangle,  and  eight-node  quadrilateral)  which  are 
available  for  use  in  the  SAMS0N2  code  (1).  SAMS0N2  is  a  dynamic 
nonlinear  two-dimensional  finite  element  computer  code  which  was 
developed  at  the  Illinois  Institute  of  Technology  Research  Institute 
(IITRI)  for  the  Air  Force  Weapons  Laboratory  (AFWL)  located  at  Kirtland 
AFB  in  Albuquerque,  New  Mexico.  SAMS0N2  is  used  to  solve 
structure-media  interface  (SMI)  problems  and  problems  which  involve 
large  displacements,  large  strains  and  nonlinear  material  behaviors. 

In  analyses  performed  previously  at  the  AFWL,  it  was  found  that  the 
solutions  obtained  by  using  the  higher-order  elements  (HOE)  were 
anomalous.  It  was  concluded  that  these  anomalous  results  were  caused  by 
inconsistencies  in  the  HOE  formulations.  Hence,  use  of  the  finite 
element  library  of  elements  for  SAMS0N2  by  the  AFWL  is  limited. 
Therefore,  the  overall  objective  of  this  research  was  to  rectify  this 
situation  by  modifying  and  correcting  as  necessary  the  HOE  formulations 
so  that  consistent  and  reliable  results  can  be  obtained.  In  addition, 
by  using  the  HOE,  it  was  believed  that  solutions  may  be  performed  more 
efficiently  with  improved  accuracy  and  fewer  elements  compared  to 
solutions  with  the  four-node  quadrilateral  continuum  element  now  in  use 


Three  tasks  were  undertaken  in  order  to  attain  the  desired 
objective  of  this  research.  The  first  task  was  to  study  previous  work 
performed  with  the  use  of  the  HOE.  The  reason  for  this  study  was  to 
acquire  background  information  as  to  where  possible  problems  may  exist 
in  the  HOE  formulations.  The  second  task  was  to  perform  finite  element 
analyses  on  various  problems  using  the  HOE  available  in  the  SAMS0N2 
code.  The  solutions  of  these  problems  were  studied  and  compared  with 
analytical  solutions  in  order  to  determine  if  the  correct  solutions  were 
obtained.  The  complexities  and  the  types  of  the  problems  which  were 
selected  were  varied  in  order  to  fully  test  the  HOE  algorithms.  The 
third  task  was  to  investigate  the  developments  of  the  HOE  algorithms  in 
the  code.  These  formulations  were  checked  for  validity  and,  if  errors 
were  identified,  corrections  for  these  errors  were  postulated.  The 
eight-node  quadrilateral  formulation  was  the  only  one  investigated 
thoroughly  since  the  eight-node  quadrilateral  is  the  most  Important  HOE 
to  the  AFWL  for  analysis  purposes.  The  five-  and  six-node  triangular 
elements  are  essentially  used  only  for  transitions  from  eight-node 
quadrilaterals  to  four-node  quadrilaterals  and  were  therefore  considered 
less  important  for  this  study. 

After  the  completion  of  these  three  tasks,  the  results  were 
compiled  and  are  discussed  in  this  report.  In  addition,  modifications 
to  the  HOE  formulation  are  suggested  and  detailed  along  with  conclusions 
based  on  these  results. 
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CHAPTER  2 

Previous  Investigations  Associated  With  the  Higher-Order  Elements 


The  only  continuum  element  currently  being  used  by  the  AFWL  is  the 
four-node  quadrilateral  (4NQ)  because  it  is  known  to  produce  correct 
solutions  and  because  AFWL  personnel  are  somewhat  concerned  about  the 
accuracy  of  the  solutions  produced  by  the  HOE.  Due  to  this  concern,  a 
few  investigations  pertaining  to  the  HOE  have  been  undertaken. 

The  content  of  this  chapter  provides  a  discussion  of  two  of  the 
prior  examinations  of  the  HOE  that  were  reported  to  the  author.  These 
investigations  were  studied  in  order  to  gain  some  insight  into  possible 
problems  associated  with  the  HOE.  Furthermore,  the  information  obtained 
from  the  studies  was  to  be  used  as  a  starting  point  for  the  current 
investigation. 

The  majority  of  the  previous  efforts  toward  investigating  the  HOE 
has  been  performed  by  Dr.  Howard  L.  Schreyer  (2)  at  the  New  Mexico 
Engineering  Research  Institute  (NMERI)  in  Albuquerque,  New  Mexico. 
Based  on  experiences  Involving  the  five-node  and  six-node  triangular 
(5NT,  6NT)  continuum  elements  and  the  eight-node  quadrilateral  (8NQ) 
continuum  element,  Schreyer  concluded  that  use  of  these  elements  has 
yielded  anomalous  results.  In  addition,  the  matrix  that  relates  stress 
to  strain  (the  B-matrix)  and  the  internal  force  vector  have  been 
investigated  at  NMERI  and  at  the  AFWL  and  have  been  accepted  as  being 
correct.  Based  on  th.s  information,  Schreyer' s  first  approach  in  his 
investigation  was  to  analyze  the  one-dimensional  wave  propagation 
problem  shown  in  Appendix  A  (Case  Al)  of  the  SAMS0N2  Users  Manual  (3) 
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using  the  8NQ  element  available  in  the  SAMS0N2  code.  By  performing  this 
analysis  he  found  that  at  the  second  time  step  in  the  solution  the 
right-side  nodes  of  the  first  element  did  not  move  as  predicted. 
Instead  he  found  the  corner  nodes  moved  to  the  left  (in  the  opposite 
direction  of  the  applied  force)  while  the  midside  node  moved  to  the 
right.  He  also  found  this  displacement  pattern  did  not  correct  with 
time.  Based  on  these  findings  he  concluded  that  the  shape  functions 
installed  in  SAMS0N2  for  the  8NQ  may  not  be  appropriate  and  that  strains 
might  be  produced  by  rigid  body  motions.  He  then  checked  the  8NQ  shape 
functions  for  rigid  body  displacement  and  rotation  and  found  no  strains 
were  produced.  After  this  check  and  with  the  results  from  the 
one-dimensional  wave  propagation  problem,  he  concluded  there  was  a 
different  reason  for  the  anomaly  and  began  reviewing  the  internal  force 
computations  and  a  modified  system  of  shape  functions.  This  last  review 
was  never  completed. 

The  only  other  reported  work  was  conducted  by  Rod  Galloway,  an 
employee  at  the  AFWt.  Doug  Seemann,  the  AFWL  project  representative  for 
the  current  Investigation  of  the  HOE,  stated  that  Rod  had  inspected  the 
HOE  formulation  In  SAMS0N2  and  found  no  obvious  errors.  Unfortunately, 
Rod's  work  was  not  documented. 

Doug  Seemann  also  expressed  some  possible  reasons  why  the  HOE 
formulation  was  not  functioning  properly  based  on  his  review  of 
Schreyer's  work  (4).  He  stated  that  tne  element  In  question  might  have 
been  unstable  due  to  the  fact  that  the  finite  element  discretization  and 
problem  input  were  poorly  defined.  He  also  believed  that  approximations 
used  in  the  element  formulation  might  be  causing  an  instability.  If  the 


problem  were  instability,  he  suggested  three  possible  solutions:  (1'  a 
more  accurate  formulation  (assuming  approximations  a^e  causing  the 
error),  (2)  more  rigid  specifications  on  how  to  use  t*e  elements  and 
(3)  a  better  scheme  for  determining  a  stable  time  step. 
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CHAPTER  3 

Finite  Element  Analyses 

For  this  study,  various  problems  were  selected,  modeled  and 
analyzed  using  the  higher-order  finite  elements  in  the  SAMS0N2  computer 
code  (1).  Most  of  the  problems  chosen  had  known  analytical  solutions  so 
that  the  results  generated  by  SAMS0N2  could  be  examined  and  evaluated 
for  errors. 

There  were  two  purposes  for  analyzing  these  problems.  First,  the 
capabilities  of  the  HOE  were  to  be  thoroughly  explored.  Second,  these 
analyses  were  to  be  used  to  discover  possible  errors  in  the  HOE 
formulation  used  in  the  SAMS0N2  code.  These  two  goals  were  obtained  by 
varying  the  complexities  and  the  types  of  problems  which  were  analyzed. 

The  8NQ  continuum  element  was  given  the  main  emphasis  in  the 
analyses  with  only  minimal  use  of  the  5NT  and  6NT  continuum  elements. 
The  4NQ  continuum  element  was  also  used  for  comparison  purposes, 
especially  for  problems  without  analytical  solutions. 

The  problem  configuration,  loading,  and  material  parameters  are 
provided  in  the  succeeding  sections  along  with  a  description  of  each 
test  problem.  Pertinent  results  are  also  provided.  Examples  of  input 
data  are  provided  in  Appendix  A  for  each  problem. 

3.1  One-Dimensional  Wave  Propagation 

A  uniform  bar  subjected  to  a  dynamic  axial  load  was  chosen  to  test 
the  solution  involving  one-dimensional  wave  propagation  which  was 
obtained  with  the  use  of  the  8NQ.  The  bar  problem  was  selected  for 
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analysis  because  Or.  Schreyer  (2)  had  used  it  in  his  previous 
investigation  of  the  higher-order  elements,  and,  therefore,  the  results 
and  observations  from  his  work  could  be  verified.  The  one-dimensional 
problem  was  also  used  in  Appendix  A  of  the  SAMS0N2  Users  Manual  (3)  to 
test  the  4NQ. 

The  geometric  configuration  of  the  uniform  bar  can  be  seen  in 

Figure  3.1  along  with  some  of  the  input  parameters  and  the  load  used  in 

the  analyses.  The  load  which  was  used  was  a  displacement  function 
applied  to  the  left  boundary  of  the  bar.  The  nodes  on  the  left  and 
right  boundaries  of  the  bar  were  fixed  in  the  vertical  (y)  direction  and 
free  in  the  horizontal  (x)  direction.  The  remainder  of  the  input  used 
for  the  analyses  is  given  in  Appendix  A.  This  input  is  identical  to 

that  in  Reference  3  for  the  4NQ.  The  input  required  for  the  8NQ  for 

this  problem  is  slightly  different  than  the  Input  for  the  4NQ.  The 
input  for  the  8NQ  solution  required  an  order  of  integration  which  was 
chosen  as  2.0  for  this  analysis.  The  bar  problem  was  analyzed 
dynamically  (undamped)  using  an  elastic  plane  strain  material  model. 

The  discretizations  used  in  the  analyses  of  the  one-dimensional 
wave  problem  were  16-4NQ  elements  with  27  nodes  and  16-8NQ  with  69  nodes 
both  divided  into  2  rows  of  8  elements.  The  4NQ  discretization  was  the 
same  as  Case  A1  mentioned  previously.  Hence,  the  same  number  of  8NQ 
elements  were  used  since  two  rows  of  elements  were  desired  with  each 
element  having  an  aspect  ratio  (ratio  of  length  of  element  to  height)  of 
1.0.  The  same  discretization  was  used  for  each  ana"  sis  since  the 
relative  improvement  in  the  results  by  increasing  or  decreasing  the 
number  of  elements  was  not  important  in  this  study. 


(fw)  iNimavidSia 


a)  Geometric  Configuration 


At  =  0.25  sec  (time  step) 

E  «  1.0  dyne/cm  (modulus  of  elasticity) 

0  A 

0  *  1.0  dyne-sec  /cm  (mass  density) 

v  *  0.0  (Poisson's  ratio) 

u  *  0.0  (damping  ratio) 

b)  Input  Parameters 


time  csecj 

c)  Applied  Load 


Figure  3.1  One-Dimensional  Wave  Propagation  Problem 
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The  results  generated  by  SAMS0N2  for  this  wave  propagation  problem 
were  compared  with  the  results  of  the  corresponding  analytical  solution 
which  can  be  found  in  either  Appendix  A  of  the  SAMS0N2  Users  Manual 
(with  u>  =  ■££)  or  Clough  and  Penzien  (5).  The  first  observation  made 
when  comparing  the  8NQ  solution  with  the  analytical  solution  was  that 
Schreyer's  statements  were  only  partially  true.  The  right-side  corner 
nodes  of  the  first  element  did  displace  to  the  left  (negative  sense) 
which  is  opposite  to  the  applied  force.  But,  contrary  to  what  Schreyer 
observed,  these  nodal  displacements  did  change  to  the  predicted 
direction  (in  the  direction  of  the  applied  force)  in  time  (15  time 
steps).  A  similar  displacement  pattern  also  occurred  for  other  nodes. 
The  magnitudes  of  the  initial  negative  displacements  had  minimal  effects 
on  the  final  results  as  can  be  seen  in  Figure  3.2.  It  has  been 
concluded  that  these  initial  negative  displacements  were  caused  by 
numerical  dispersion  and  were  inherent  in  the  8NQ  shape  functions  and 
element  formulation. 

Additional  results  can  be  seen  in  Figures  3. 3-3. 7.  The  results  for 
the  4NQ  and  the  analytical  solution  are  identical  to  those  shown  in 
Appenaix  A  of  the  SAMS0N2  Users  Manual.  The  8NQ  results  compare  with 
the  analytical  solution  quite  well,  except  for  the  oscillatory  numerical 
dispersion  present.  Only  one  displacement  graph  is  given  which  is 
representative  of  all  the  other  displacement  graphs.  All  displacement 
graphs  showed  similar  degrees  of  correlation  between  the  8NQ  and  the 
analytical  solution  results.  The  displacement  and  stress  graphs  have 
the  best  agreement  between  the  8N0  solution  and  the  analytical  solution, 
while  the  velocity  graphs  displayed  a  greater  amount  of  dispersion. 
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Flqure  3.7  One-Dimensional  Have  Propagation  Problem:  Stress  vs. 


One  conclusion,  based  on  the  results  generatec  by  SAMSON2,  was  that 
the  8NC  element  is  satisfactory  for  use  in  the  one-dimensional  wave 
prooagation  problem  for  which  a  dynamic  elastic  solution  was  performed. 
The  only  discrepancies  in  any  of  the  results  which  should  be  noted 
existed  in  the  values  for  the  deflections,  velocities,  and  stresses  in 
the  y-direction  at  points  where  the  4NQ  solution  and  the  8NQ  were 
compared.  However,  the  magnitudes  of  these  discrepancies  were 
insignificant. 

Some  other  conclusions  that  were  obtained  include:  1)  the 

formulation  for  the  8NQ  has  some  inherent  problems,  and  2)  the  solution 
obtained  with  the  4NQ  was  more  efficient  than  that  obtained  with  the 
8NQ.  The  CPU  time  for  the  4NQ  solution  was  7.20  secs,  compared  to  24.33 
secs,  for  the  8NQ. 


3.2  Cantilever  Beam 

A  problem  involving  a  cantilever  beam  was  solved  in  order  to  verify 
the  formulation  for  the  4NQ  which  appears  in  both  the  SAMS0N2  Users 
Manual  and  the  STEALTH  and  SAMS0N2  Verifications  Manual  (6).  Hence, 
this  beam  problem  is  used  to  test  the  performance  of  the  5NT,  6NT,  and 
8NQ  elements  in  flexure  when  analyzed  both  statically  and  dynamically. 

The  geometric  configuration  of  the  cantilever  beam  can  be  seen  in 
Figure  3.8  along  with  a  sample  of  the  pertinent  input  data  and  loads 
used  in  the  analyses.  The  nodes  on  the  left  boundary  were  fixed  in  both 
the  x-  and  y-directions.  The  maximum  stable  time  step  for  the  uNQ 
discretization  was  obtained  by  trial  and  error  (based  on  successive 
unsuccessful  computer  executions)  until  a  satisfactory  solution  was 
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obtained.  Once  chosen,  the  same  time  step  was  used  for  all  of  the 
analyses.  Mass  proportional  damping  (Cj)  was  used  to  obtain  the  static 
solutions  by  dynamic  relaxation  techniques  (3).  The  beam  was  subjected 
to  a  static  and  a  dynamic  uniformly  distributed  pressure  load.  A  short 
rise  time  was  used  for  the  static  analyses  in  order  to  reduce  the  time 
required  for  convergence  to  the  solution.  The  material  law  for  a  state 
of  biaxial  elastic-plastic  plane  stress  was  used  in  the  analyses  with 
Poisson's  ratio  equal  to  0.0.  An  order  of  integration  of  2.0  was  chosen 
for  the  5NT,  6NT,  and  8NQ.  Examples  of  the  remainder  of  the  input  data 
used  for  the  static  analyses  are  contained  in  Appendix  A.  Only  two 
changes  were  made  to  this  input  data  in  order  to  execute  the  dynamic 
analyses.  The  first  change  was  the  removal  of  the  mass  proportional 
damping  factor  from  the  input.  The  second  change  was  to  eliminate  the 
short  rise  time  in  the  static  loading  curve  in  order  to  make  it  an 
instantaneous  dynamic  load. 

The  finite  element  discretizations  used  for  the  analyses  of  the 
cantilever  beam  are  shown  In  Figure  3.9a,  b,  c  and  d.  The  4NQ  grid  was 
based  on  the  example  problems  in  the  SAMS0N2  Users  Manual.  The  element 
dimensions  were  selected  to  be  1  ft  x  1  ft  in  order  to  maintain  an 
aspect  ratio  equal  to  one  as  recommended  in  Reference  6.  The  8NQ 
discretization  was  selected  such  that  one  8NQ  element  replaced  four  4NQ 
elements.  This  ratio  of  one  8NQ  to  four  4NQs  was  chosen  so  that  the 
same  number  of  stress  and  strain  calculations  for  the  two  discretization 
would  be  performed.  Hence,  the  efficiencies  of  the  two  solution  schemes 
could  be  compared.  The  aspect  ratio  for  the  8NQ  was  also  1.0.  The 
discretizations  for  both  the  5NT  and  6NT  contained  32  elements. 


c)  5NT  Elements 
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d)  6 NT  Elements 


Cantilever  Beam  Problem:  Finite  Element  Discretizations 
(Continued) . 
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Therefore,  the  number  of  5NT  and  6NT  elements  used  in  the  analyses  was 
twice  the  number  of  8NQ  elements  and  one-half  the  number  of  4NQ 
elements.  Fewer  triangular  elements  ( 5NT  and  6NT)  were  selected 
compared  to  the  number  of  4NQ  because  the  order  of  integration  was 
greater  for  the  5NT  and  6NT  than  for  the  4NQ.  Once  selected,  the 
discretizations  used  in  the  analyses  were  never  changed. 

The  results  of  the  static  analyses  for  each  of  the  four  element 
discretizations  were  compared  to  the  corresponding  analytical  solutions. 
These  solutions  were  obtained  with  the  use  of  basic  principles  of 
mechanics  (7).  The  maximum  y-displacements  (y_,v)  obtained  from  the 
computer  analyses  were  compared  to  the  theoretical  value  computed  from 


the  following  equation: 


wl4  3wL2 

W  *  -  srr  ‘tft 


(3.1) 


where. 


w  =  200.0  lb/ft 
L  =  16.0  ft 
E  =  4176.0  x  106  psf 
1  =  10.667  ft4 
A  =  8.0  ft2 
therefore, 

ymax  3  -°-3908  x  10’4  ft  • 

This  value  for  the  maximum  displacement  was  computed  at  x  =  L  and 
included  shear  deformation  effects.  Table  3.1  shows  the  comparison 
between  the  results  generated  by  SAMS0N2  for  the  maximum  y-displacement 


and  the  value  computed  using  the  analytical  solution. 
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Table  3.1 


Comparison  Between  the  Results  Generated  by  SAMSON 2  for 
the  Maximum  Y-Displacement  and  the  Corresponding 
Analytical  Solution. 


Element 

Type 

SAMS0N2  . 

y_,  (ft  x  10"  ) 
max 

Analytical 
Solution  . 

Wft  *  10  > 

Percent  Difference 
Based  On 

Analytical  Solution 

4NQ 

-0.4067  (Node  51) 

-0.3908 

4.07% 

8NQ 

-0.3861  (Node  43) 

-0.3908 

1 . 20% 

5NT 

-0.3535  (Node  35) 

-0.3908 

9.54% 

6NT 

-0.3858  (Node  51) 

-0.3908 

** 

CO 

CNJ 

The  results  for  the  4NQ,  8NQ,  and  6NT  elements  as  seen  from  Table  3.1 
differed  only  slightly  when  compared  to  the  value  obtained  from  the 
analytical  solution.  The  5NT  element  results  were  not  acceptable  as 
they  varied  from  the  analytical  solution  by  more  than  10%.  A  comparison 
between  the  analytical  and  SAMS0N2  solutions  for  normal  stresses  (o  ) 
was  also  performed.  The  normal  stress  values  used  for  comparison  were 
computed  analytically  using  the  following  equation: 

My  (100x2  -  3200x  +  25,600)y 


moment  in  the  beam  at  any  distance  x  (measured  from  the  left 
end  of  beam)  0.0  _<  x  <  16.0, 

distance  from  the  neutral  axis  (N.A.)  of  the  beam 
-2.0  <  y  <  2.0,  and 

moment  of  inertia  of  the  cross-sectional  area  of  the  beam  with 
respect  to  the  x-axis. 
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Table  3.2  shows  some  of  the  comparisons  made  between  the  computed  ncrma1 
stresses  from  Equation  3.2  to  those  computed  using  SAMS0N2.  As  before, 
the  results  for  the  4NQ,  8NQ,  and  6NT  elements  compared  well,  but 

results  for  the  5NT  did  not.  Normal  stresses  for  other  elements  in  the 
meshes  were  compared,  but  are  not  shown  here.  Similar  observations  for 
the  stresses  in  these  other  elements  were  made  witn  regard  to  the  four- 
different  element  discretization  solutions. 

The  results  from  the  dynamic  analyses  were  also  compared  to  the 
corresponding  analytical  solution  and  are  shown  in  Figures  3.10  to  3.17 
and  Tables  3.3  and  3.4.  The  analytical  solution  for  the  displacement 
response  at  the  free  end  of  the  cantilever  beam  was  generated  using 

Equations  3.3,  3.4,  and  3.5  (8).  These  equations  are  for  the  elementary 
case  and  do  not  include  shear  distortion  and  rotary  inertia  or 

axial-force  effects.  These  three  quantities  were  omitted  to  simplify 

the  solution. 

V(x,t)  =  r  «n(x)  Y(t)  (3.3) 

n-1 

where, 

V(x,t)  =  displacement  response  at  any  x  for  any  time  t, 

$n(x)  =  vibration  shape  at  any  x  for  mode  n,  and 

Yn(t)  *  amplitude  of  vibration  for  mode  n  at  any  time  t. 

$n(x)  =  (cosh  anx  -  cos  apx)  -  Cn(sinh  apx  -  sin  anx)  (3.4) 


where , 


Table  3.2  A  Comparison  Between  the  SAMS0N2  Results  and  the  Values  Obtained  by  Use  of  Equation  3.2 
for  Normal  Stresses  (o  ). 
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x  =  distance  from  fixed  end,  and 

a:  =  1.875/L  ,  a2  =  4.694/L  ,  and  a3  =  7.855/L. 

Yn(t}  =  ~r  (zrr^1  -  cos  wnt) 
m  n 

where, 


(3.5) 


PQ  =  constant  value  for  the  uniformly  distributed  load  (lb/ft), 
m  =  constant  value  for  mass  per  unit  length  (lb-sec  /ft  ), 


2 

un 


t  = 


L 


$  (x)dx 

—  -  ;  I.  =  0.7830  ,  I,  =  0.4340  ,  I, 

*£(x)dX  1  2  3 

El 

—  =  frequency  of  vibration  for  mode  n; 


0.2589, 


m 

Wj  *  256.18  ,  u)£  *  1605,59  ,  a  4,496.16,  and 
time  in  sec. 
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The  first  three  modes  of  vibration  (n=l,2,3)  were  used  to  predict  the 
displacement  responses  of  the  cantilever  beam  as  shown  in  Figures  3.10 
and  3.11.  As  shown  in  these  figures,  the  analyses  using  the  4NQ,  8NQ, 
and  6NT  elements  all  overestimate  the  results  from  the  analytical 
solution.  In  addition,  the  response  of  these  elements  has  a  longer 
period  of  vibration  relative  to  that  of  the  exact  solution.  The 
opposite  is  true  of  the  5NT  element  response,  i.e.,  the  5NT  solution 
underestimates  displacements  and  has  a  shorter  period  of  vibration.  The 
4NQ  and  the  exact  solution  correspond  to  what  is  presented  in 
Reference  3.  A  numerical  comparison  of  the  average  peak  values  for  the 
displacements  for  each  curve  in  each  figure  can  be  seen  in  Table  3.3 
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Table  3.3  Comparison  of  Average  Values  of  Peak  Displacement  and 
Average  Periods  of  Oscillation  for  the  Curves  Shown  in 
Figures  3.10  and  3.11. 


Solution 

Type 

Average  Value 
of  Peak 
Displacement 

(ft  x  10‘4) 

1 

Percent 
Difference 
Based  on 
Analytical 
Solution 

2 

Percent 

Difference 

When  Compared  to 

2  Times  Static 
Displacement 

(0.7816  x  10"4) 

3 

Average 
Period  of 
Oscillation 
(sec) 

Analytical 

-0.7403 

— 

— 

0.0248 

4NQ 

-0.8238 

11.3* 

5.4% 

0.0259 

8NQ 

-0.7793 

5.3% 

0.3% 

0.0255 

5NT 

-0.7134 

3.6% 

8.7% 

0.0219 

6N.T 

-0.7810 

5.5% 

0.1% 

0.0254 

along  with  the  average  periods  of  oscillation.  In  this  table,  the  first 
comparison  was  made  between  the  average  peak  value  for  the  exact  curve 
and  each  of  the  other  four  curves.  These  results  were  different  than 
those  obtained  from  the  static  analyses  in  that  the  5NT  element  solution 
exhibited  the  best  correlation  to  the  analytical  solution.  It  was 
concluded  that  this  good  correlation  resulted  from  neglecting  shear 
distortion  in  the  analytical  solution.  Another  comparison  was  performed 
between  a  dynamic  displacement  response  value  equal  to  two  times  the 
static  displacement  and  the  average  peak  values  from  each  of  the  other 
four  curves.  The  conclusions  from  this  comparison  were  similar  to  those 
from  the  static  solution,  and,  hence,  the  results  from  the  4NQ,  8NQ,  and 
6NT  element  solutions  are  acceptable.  The  periods  of  vibration  for 
these  three  element  types  were  also  in  good  agreement  with  those  of  the 
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exact  solution.  The  5NT  element  solution  had  the  largest  discrepancies 
in  columns  2  and  3  of  Table  3.3,  but  the  solution  was  acceptable  based 
on  column  1  in  which  shear  distortion  was  neglected.  The  dynamic  normal 
stress  (o  )  responses  generated  for  the  four  element  types  were  also 
studied.  The  following  equation  was  used  to  compute  the  stress  response 
values  that  were  used  to  compare  with  the  output  from  the  analyses  using 
the  four  different  element  types: 


®  (x.t)  -  Ey  I  *"(x)  Y  (t) 
x  n=l  n 


(3.6) 


where, 


•;(x)  =  a^  [(cosh  anx  +  cos  anx)  -  Cn(sinh  anx  +  sin  anx ) ] . 

Curves  showing  the  variation  in  the  maximum  stress  and  the  stress  at 
midspan  (same  locations  and  elements  as  in  the  static  analyses)  for  each 
element  type  are  given  in  Figures  3.12  to  3.17.  The  results  are  similar 
to  those  from  the  static  analyses  and  the  displacement  response 
comparisons.  The  4NQ,  8NQ,  and  6NT  element  solutions  agree  well  with 
the  exact  solution,  but  the  5NT  solution  does  not  agree  with  the 
analytical  solution  as  well.  Table  3.4  shows  the  comparison  between  the 
peak  stress  values  for  each  element  type  relative  to  the  exact  solution. 
The  6 NT  and  8NQ  higher-order  elements  agree  with  the  exact  solution  the 
best.  The  4NQ  element  solution  is  acceptable,  but  the  5NT  element  gives 
the  worst  correlation  to  the  exact  solution. 
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The  execution  times  (CPU  times)  were  recorded  for  each  analysis  of 
the  cantilever  beam  in  addition  to  the  results  presented  previously. 
These  CPU  times  were  used  to  determine  the  relative  efficiencies  of  the 
solutions  using  the  four  element  types.  Table  3.5  shows  a 
representative  sample  of  the  CPU  times  obtained  from  the  analyses.  As 
shown,  the  4NQ  element  solution  was  two  to  three  times  more  efficient 
than  the  three  HOE  solutions.  These  significant  increases  in  CPU  time 
for  the  HOE  solutions  as  compared  to  the  4NQ  solution  were  attributed  to 
the  following:  (1)  more  terms  are  included  for  each  calculation 
involving  the  shape  functions  for  the  HOE  due  to  the  increased  number  of 
nodes  per  element,  (2)  more  subroutines  are  used  in  the  HOE  solution 
schemes,  and  (3)  the  5NT  and  6NT  element  solutions  involve  twice  as  many 
stress  calculations  as  the  4NQ  element  solution.  It  should  be  noted 
that  these  three  conditions  more  than  offset  the  time  saved  by  the 
reduction  in  the  number  of  the  nodes  resulting  from  the  use  of  the  HOE. 
Based  on  the  third  condition,  if  the  same  number  of  stress  calculations 
were  performed,  an  increase  in  CPU  time  of  between  50  and  75  percent 
would  be  expected  for  the  5NT  and  6NT  element  solutions  when  compared  to 
the  4NQ  solution. 

Some' general  conclusions  which  are  based  on  the  results  presented 
in  the  previous  paragraphs,  are  as  follows.  The  6NT  and  8NQ 
higher-order  elements  and  the  4NQ  element  give  correct  solutions  for 
both  static  and  dynamic  analyses  for  a  problem  involving  flexural 
response  with  an  elastic  plane  stress  material  law.  Solutions  involving 
these  two  HOE  showed  more  accurate  results  with  fewer  nodes  and  elements 
than  the  4NQ  element  solutions,  but  were  much  less  efficient  (more 
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Table  3.5 

CPU  Time  Comparisons  for  the 
the  Analyses  of  the  Cantilever 

Four  Element  Types  Used  in 
Beam. 

Element  Type 

CPU  Time 

Percent  Difference  Based 

(secs) 

on  4NQ  Solution 

4NQ 

105.0 

— 

8NQ 

224.1 

113.4% 

5NT 

317.0 

201.9% 

6NT 

339.8 

223.6% 

expensive).  The  decrease  in  the  accuracy  of  the  4NQ  element  solution  is 
due  to  the  under-integration  inherent  in  the  formulation  of  this  element 
as  it  exists  in  the  SAMS0N2  code  (3).  The  5NT  element  results  varied 
greatly  with  respect  to  results  from  the  other  three  elements  for  both 
the  static  and  dynamic  cases,  except  for  the  maximum  dynamic  displace¬ 
ment  which  compared  well  to  results  from  an  analytical  solution  which 
neglects  the  effects  of  shear  deformation.  These  poor  results  were 
believed  to  be  due  to  the  fact  that  the  5NT  element  is  not  an  entirely 
linear  strain  triangle  because  it  only  has  two  quadratic  displacement 
sides.  This  fact  could  cause  the  5NT  element  to  behave  in  a  manner 
similar  to  a  constant  strain  triangle  in  that  a  soft  response  might 
occur  which  could  only  be  corrected  by  increasing  the  number  of  elements 
used  in  the  analysis  (3).  Another  possible  cause  of  the  poor  results 
exhibited  by  the  5NT  element  was  the  position  of  the  midside  nodes  in 
the  5NT  elements.  For  the  analyses  which  were  performed,  the  midside 
nodes  were  positioned  on  the  hypotenuse  and  on  one  side  of  the  triangle. 
A  different  discretization  was  devised  with  the  two  midside  nodes  being 


on  the  two  sides  of  the  triangle  and  not  on  the  hypotenuse.  A  solution 
using  this  different  discretization  involving  the  nodes  on  the  two  sides 
was  attempted,  but  no  output  was  obtained  because  of  an  incompatibility 
(zero  area  was  computed  by  the  program)  between  the  input  and  the 
corresponding  formulation  of  the  5NT  element  which  could  not  be 
resolved.  A  final  note  on  the  5NQ  element  is  that  this  element  did 
produce  some  acceptable  results  when  compared  to  the  analytical 
solution.  Therefore,  it  might  be  possible  to  make  use  of  this  element 
for  the  transition  from  an  8NQ  element  to  a  4NQ  element  as  it  was 
intended. 

3.3  Soil -Structure  Interaction 

Based  on  the  results  for  the  HOE  solutions  from  the  previous  two 
sections,  especially  the  results  for  the  8NQ  element,  a  more  complex 
problem  was  chosen  in  order  to  further  check  the  performance  of  the  8NQ 
element.  A  problem  involving  soil-structure  interaction,  the  type  of 
problem  for  which  SAMS0N2  was  designed,  was  chosen.  For  this  problem, 
no  analytical  solution  Is  available  and  all  of  the  8NQ  results  were 
compared  to  those  obtained  by  using  the  4NQ.  Although  this  comparison 
may  not  appear  to  be  very  practical,  it  gave  a  qualitative  measure  of 
the  performance  of  the  8NQ  element  for  this  problem.  The  4NQ  has  been 
shown  by  personnel  at  the  AFWL  to  be  reliable  and  accurate  for  similar 
soi 1-structure  interaction  problems  when  compared  to  actual  test  data. 

The  configuration  of  the  soil-structure  problem  is  shown  in 
Figure  3.18a,  b  and  c  along  with  some  input  data  and  the  primary  load 


a)  Problem  Configuration 


1.0E-6  sec 

. 2516E-3  lb-sec2/in4 
. 173E-3  lb-sec2/1n4 
. 42E7  psi 

. 158E7  psi 


(time  step) 

(mass  density  of  concrete) 

(mass  density  of  soil) 

(modulus  of  elasticity  of 
concrete) 

(modulus  of  elasticity  of  soil) 


0.24 

0.38 

0.01 


(Poisson's  ratio  for  concrete) 
(Poisson's  ratio  for  soil) 
(damping  ratio) 
b)  Input  Parameters 


Figure  3.18  Soil -Structure  Interaction  Problem. 


2 


(nuclear  airblast)  used  in  the  analyses.  This  problem  is  a  scaled 
version  of  one  which  was  analyzed  and  experimentally  tested  by  the  AFWL. 
The  structure  is  a  fiber  reinforced  concrete  cylinder  surrounded  by  Yuma 
soil  separated  by  a  structure-media  interface.  The  input  data 
pertaining  to  the  two  materials  were  obtained  from  AFWL  material  models. 
Additional  input  data  are  contained  in  Appendix  A.  Also,  the  slideline 
data  are  from  the  AFWL  bilinear  failure  surface  model  for  Yuma  soil  and 
fiber  reinforced  concrete.  An  integration  order  of  2.0  was  used  in  the 
8NQ  analyses.  The  load  curve  was  obtained  from  an  actual  test  performed 
by  the  AFWL  and  was  estimated  for  use  in  the  analyses  with  the  use  of 
the  NMERI  Speicher/Brode  Nuclear  Airblast  Curve  algorithm  based  on  a 
peak  pressure  of  50  ksi  and  a  yield  of  225  kt. 

In  the  analysis  of  the  soil -structure  problem,  the  nodes  on  the 
bottom  boundary  of  the  soil  field  were  not  allowed  to  displace  in  the 
y-direction,  while  those  on  the  right  boundary  could  not  displace  in  the 
x-di recti  on.  The  lower  right  corner  node  was  fixed  in  both  the  x-  and 
y-directions.  The  problem  was  analyzed  dynamically  using  an 
axi symmetric  material  law  developed  by  the  AFWL.  The  law  is  called  the 
AFWL  "engineering"  model  and  was  primarily  designed  to  model  soil 
behavior.  It  is  defined  by  a  piecewise  linear  hydrostat  (a  plot  of 
hydrostatic  pressure  vs.  volumetric  strain)  for  which  the  critical 
identifying  points  are  given  as  part  of  the  input  data  (3).  Strain 
softening  and  dilation  cannot  be  modeled  by  the  AFWL  "engineering"  law. 
Therefore,  the  hydrostat  curve  must  be  monotonlcally  increasing  (strain 
hardening)  (6).  It  is  important  to  note  that  the  soil-structure  problem 
should  be  analyzed  in  three  dimensions  in  order  to  obtain  the  most 
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accurate  results.  However,  SAMS0N2  is  a  two-dimensional  code,  and, 
therefore,  this  problem  was  analyzed  axisymmetrically  with  SAMS0N2. 

The  soil-structure  interaction  problem  was  discretized  separately 
into  4NQ  and  8NQ  elements  for  the  analyses.  The  initial  finite  element 
mesh  for  the  4NQ  consisted  of  216  -  6  in.  x  6  in.  elements  (204  soil 
elements,  12  concrete  elements)  and  the  initial  mesh  for  the  8NQ 
consisted  of  54  -  12  in.  x  12  in.  elements  (51  soil  elements,  3  concrete 
elements).  These  meshes  contained  relatively  large  elements  (a  coarse 
discretization)  in  order  to  keep  the  CPU  times  as  short  as  possible, 
especially  for  the  8NQ  solutions,  and  because  only  qualitative  results 
were  desired  for  comparison  purposes.  The  8NQ  discretization  was  chosen 
such  that  one  8NQ  element  replaced  four  4NQ  elements.  Therefore,  a  CPU 
time  comparison  could  be  made  as  was  done  for  the  cantilever  beam 
problem. 

A  qualitative  comparison  was  performed  between  the  two  solutions. 
In  making  this  comparison,  the  principal  values  considered  were  the 
displacements  and  accelerations  in  the  y-direction  for  all  the  nodes  in 
the  concrete  cylinder  as  well  as  for  a  representative  number  of  nodes  in 
the  soil.  Also,  the  stresses  in  the  y-direction  ( o  )  were  considered 
important  for  all  elements  in  the  concrete  and  for  a  typical  number  of 
soil  elements.  The  values  for  these  three  output  parameters  were 
studied,  both  qualitatively  with  regard  to  their  variation  in  response 
with  time  and  numerically. 

The  comparisons  for  the  initial  (coarse)  meshes  rendered  very  few 
acceptable  results.  In  general,  the  results  obtained  from  the  8N0 
element  solution  using  the  initial  coarse  mesh  in  the  concrete 
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(Figure  3.19b)  exhibited  a  significant  amount  of  variation.  For 
instance,  the  upper  right  node  in  the  structural  mesh  (node  202)  had  a 
y-displacement  value  three  times  that  of  the  other  two  adjacent  nodes 
(nodes  200  and  201)  at  the  end  of  the  tenth  time  step  (the  first  time 
step  for  which  output  was  obtained).  This  variation  was  thought  to  be 
due  to  one  of  three  conditions:  (1)  the  mesh  which  was  used  was  too 
coarse,  (2)  instability  existed  in  the  solution,  or  (3)  the  slideline 
calculations  for  the  8NQ  element  were  in  error.  The  4NQ  element 
solution,  in  contrast  to  the  8NQ  solution,  had  very  consistent  results. 
In  fact,  the  top  three  nodes  in  the  coarse  4NQ  finite  element  mesh  of 
the  concrete  cylinder  (nodes  242,  249,  and  256  in  Figure  3.19a)  had 
equal  values  for  both  the  y-displacements  and  y-accelerations  at  the  end 
of  the  tenth  time  step.  The  stresses  for  the  4NQ  solution  were  also 
consistent,  while  those  for  the  8NQ  were  erratic  due  to  the  unequal 
y-displacements  and  y-accelerations  of  the  upper  nodes.  A  comparison  of 
the  values  for  the  y-displacements,  y-accelerations,  and  stresses  in  the 
y-direction  (oy)  was  performed  at  the  end  of  the  tenth  time  step.  The 
y-displacement  and  y-acceleration  values  for  nodes  242  and  249  in  the 
4NQ  mesh  were  twice  those  of  nodes  200  and  201  in  the  8NQ  mesh,  while 
the  values  for  node  256  in  the  4NQ  mesh  were  significantly  less  than 
those  for  node  202  in  the  8NQ  mesh.  A  comparison  of  stresses  between 
the  top  element  in  the  8NQ  mesh  (element  54)  and  the  top  four  elements 
in  the  4NQ  mesh  (elements  209,  210,  215  and  216)  showed  an  acceptable 
correlation  for  only  one  stress  value.  The  stress  at  the  upper  right 
integration  point  of  element  54  was  within  12?  of  the  stress  in  element 
216.  The  normal  stress  (o  )  also  compared  well  at  this  location.  The 
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shearing  stresses  in  the  8NQ  solution  at  the  tenth  time  step  were 
significantly  higher  than  those  from  the  4NQ  solution.  These  large 
differences  in  the  shearing  stresses  as  well  as  for  the  other  stresses 
resulted  from  the  large  numerical  discrepancies  in  the  y-di splacements 
at  the  top  of  the  8NQ  mesh.  The  other  nodes  in  the  8NQ  mesh  behaved 
much  the  same  as  the  top  row  of  nodes,  but  exhibited  one  additional 
trend  which  was  previously  discovered  in  the  one-dimensional  wave 
propagation  problem.  Many  of  the  nodes  located  below  the  second  row  of 
nodes  in  the  8NQ  mesh  (nodes  198  and  199)  displaced  upward  (positive 
direction)  which  is  opposite  to  the  direction  of  the  applied  blast  load, 
while  the  corresponding  nodes  in  the  4NQ  all  displaced  downward 
(negative  direction). 

The  8NQ  solution  was  corrected  with  time  as  in  the  wave  propagation 
problem,  but  only  to  a  certain  extent.  Most  of  the  nodes  in  the 
structure  did  reverse  directions  by  the  end  of  400  time  steps.  The 
y-displacement  values  for  nodes  200,  201  and  202  did  tend  toward  the 
same  value,  but  the  y-dlsplacement  for  node  202  was  still  greater.  The 
y-di splacements  for  the  nodes  on  the  right  vertical  boundary  (slideline 
nodes  202,  199,  197,  194,  102,  189,  and  187)  were,  in  most  instances, 
greater  than  the  y-di splacements  for  the  corresponding  nodes  to  the  left 
of  the  slideline.  Stress  values  (oy)  at  similar  levels  within  the 
structural  mesh  were  approximately  equal  (consistent  at  these  levels)  at 
certain  times  during  the  solution,  but  this  consistency  was  not  observed 
on  a  regular  basis.  Oscillations  (wave  propagations)  were  apparent 
throughout  the  solution  as  was  shown  in  the  plots  of  y-accelerations , 
stresses  in  the  y-direction  (o^),  and  y-displacements  in  the  SAMS0N2 


output.  Results  from  the  soil  elements  were  consistent  for  the  entire 
time  of  solution. 


Quantitative  comparisons  of  results  between  the  4NQ  solution 
(solution  remained  consistent  with  the  nodes  being  almost  equal  for  each 
individual  row,  but  those  on  the  slideline  displaced  slightly  more  than 
the  others)  and  the  8NQ  solution  were  performed  at  every  tenth  time  step 
in  addition  to  a  qualitative  overall  comparison.  Some  numerical 
comparisons  showed  acceptable  results,  while  the  majority  of  the  results 
compared  poorly.  The  y-displacements  for  the  top  two  rows  of  nodes  in 
the  soil  were  within  81  of  each  other  at  t  *  .2E-3  (200  time  steps),  but 
deviated  by  up  to  37%  at  t  *  .4E-3  (400  time  steps),  respectively.  The 
best  soil  stress  correlations  differed  by  17%.  Stresses  in  the  concrete 
cylinder,  particularly  a^,  matched  well  at  certain  times,  but  varied 
with  the  propagations  of  the  waves  through  the  structures.  Therefore, 
this  correlation  of  stresses  might  have  been  by  coincidence  more  than 
anything  else.  The  large  shear  stress  discrepancies  discovered  at  the 
tenth  time  step  were  significantly  reduced  with  time.  Qualitatively, 
both  solutions  displayed  the  wave  propagation  phenomenon  throughout  the 
structure,  but  the  8NQ  solution  was  trailing  that  of  the  4NQ  as  if  the 
wave  speeds  in  the  two  were  quite  different.  The  plots  of  the 
y-displacements  in  the  SAMS0N2  output  for  the  4NQ  solution  showed  the 
wave  propagation  interaction  more  than  those  for  the  8NQ  solution.  In 
general,  the  nodal  displacements  for  the  corresponding  rows  are 
qualitatively  similar  such  that  the  displacements  are  larger  for  nodes 
closer  to  the  slideline.  The  values  for  the  nodal  y-displacements  from 
the  4NQ  solution  are  greater  than  those  from  the  8NQ  solution,  except  in 


the  top  row,  even  though  at  certain  times  there  was  good  agreement 
between  corresponding  values. 

Many  variations  in  the  solutions  were  attempted  in  order  to  improve 
the  results  stated  in  the  previous  paragraphs  for  the  initial  4NQ  and 
8NQ  coarse  discretizations.  The  variations  to  the  input  were  as 
fol 1 ows : 

(1)  energy  error  output  was  requested  for  the  two  solutions 
to  determine  if  unstable  conditions  were  present, 

(2)  the  slave  nodes  were  reordered  in  the  slideline  in  an 
attempt  to  eliminate  the  erratic  behavior  present  in  the 
8NQ  solution, 

(3)  reduced  loadings  were  applied  to  reduce  the  possibility 
of  instabilities  occurring  in  the  solution, 

(4)  finer  mesh  discretizations  of  the  concrete  cylinder  were 
used  in  the  4NQ  and  8NQ  solutions,  and 

(5)  the  order  of  integration  was  Increased  to  3.0  for  the  8NQ 
solution. 

The  results  from  requesting  the  energy  error  output  are  shown  in 
Figures  3.20  and  3.21.  These  figures,  if  assumed  reliable  based  on  the 
findings  of  Berglund  and  Rudeen  in  Reference  6,  show  that  both  solutions 
are  unstable  because  the  energy  error  exceeds  1%  (1,  9).  In 
Reference  6,  it  was  pointed  out  that  the  energy  terms  calculated  in 
SAMS0N2  are  either  in  error  or  Incomplete  and  need  to  be  corrected.  In 
the  current  study,  assuming  the  results  from  the  Figures  3.20  and  3.21 
to  be  reliable,  the  possibility  of  excessive  distortion  was  investigated 
in  order  to  determine  if  the  8NQ  formulation  was  still  valid  (11).  No 
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apparent  excessive  distortion  was  located  during  the  current  study. 
Upon  investigating  the  possibility  of  excessive  distortion,  a  study  of 
time  step  sizes  versus  changes  in  the  solutions  was  performed.  The  time 
step  size  was  halved,  and  then  halved  again.  No  apparent  changes 
occurred  in  either  of  these  two  solutions.  Therefore,  the  solutions 
were  considered  to  be  stable,  especially  since  Figures  3.20  and  3.21 
show  the  4NQ  solution  to  be  more  unstable  than  that  for  the  8NQ  for  the 
same  time  step  size.  It  was  found  in  the  execution  of  the  cantilever 
beam  problem  that  the  stable  time  step  for  the  8NQ  solution  was  less 
than  that  which  could  have  been  used  in  the  4NQ  solution.  An  additional 
note  to  this  information  is  that  the  curves  for  each  of  the  energy 
components  (kinetic,  internal,  and  external  work)  displayed  the  same 
patterns. 

The  second  variation  in  the  input  that  was  performed  in  an  attempt 
to  improve  the  8NQ  results  when  compared  to  those  obtained  for  the  4NQ, 
was  a  change  to  the  order  of  the  slideline  nodes.  This  change  was 
performed  due  to  a  phone  conversation  with  personnel  of  the  AFWL  (4). 
The  actual  change  was  just  a  reversal  of  the  order  of  the  input  for  the 
slave  nodes.  Instead  of  inputting  the  nodes  from  top  to  bottom,  they 
were  input  from  bottom  to  top  so  that  the  first  slave  node  was  adjacent 
to  the  last  master  node  on  the  slideline.  The  result  of  this  change  was 
very  minor  and  affected  mostly  the  lower  nodes  of  the  structural  meshes. 
No  relative  improvement  in  results  was  obtained. 

The  third  variation  was  the  application  of  a  reduced  load  in  order 
to  lessen  the  chance  of  instabilities  occurring  as  well  as  to  make  the 
solution  perform  more  elastically  with  little  or  no  plastic  flow.  The 
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first  load  variation  was  the  blast  load  reduced  by  a  factor  of  500  (peak 
pressure  of  100  psi  in  place  of  50,000  psi).  The  results  obtained  for 
the  initial  structural  meshes  with  this  applied  load  were  no  better  than 
those  for  the  originally  applied  airblast  load.  The  same  trends  in  the 
nodal  displacements  existed  in  these  solutions  as  in  the  previous  ones. 
Another  solution  was  executed  with  the  100  psi  load  reduced  to  a  peak 
value  of  1  psi  and  twice  the  run  time  used  in  each  of  the  previous 
solutions.  The  solution  time  was  doubled  in  order  to  determine  if  the 
8NQ  solution  would  converge  further.  The  results  from  the  solution 
using  the  1  psi  load  and  the  double  run  time  compared  well  with  those 
obtained  with  the  use  of  the  100  psi  load,  such  that  at  time  t  *  . 4E-3 
(the  total  run  time  for  the  100  psi  solution)  all  results  differed  by  a 
factor  of  100  which  was  the  same  ratio  between  the  two  applied  loads. 
These  results  showed  that  the  two  solutions  which  were  obtained  with  the 
use  of  the  reduced  loads  were  both  stable  and  in  the  elastic  range  as 
was  desired  in  these  analyses.  The  remainder  of  the  results  (from  t  = 
.4E-3  to  .8E-3)  for  the  1  psi  peak  load  revealed  the  following: 

(1)  The  nodes  on  the  slideline  in  the  two  solutions  no  longer 
had  the  largest  values  for  the  y-displacement,  but 
instead  the  displacements  got  larger  in  proceeding  toward 
the  left  side  of  the  structural  mesh  along  a  row.  This 
behavior  was  opposite  to  that  found  for  previous 
solutions.  The  slideline  nodes  acted  as  if  they  had 
reached  a  limiting  state.  This  new  displacement  pattern 
was  found  throughout  the  entire  4NQ  mesh  and  in  the  top 
two-thirds  of  the  8NQ  mesh. 
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(2)  Some  better  comparisons  were  obtained  between  the  two 
solutions  for  values  of  nodal  displacements  and  elemental 
stresses.  But,  as  before,  these  good  results  were  only 
obtained  for  particular  locations  and  not  throughout  the 
structural  meshes. 

(3)  The  8NQ  results  were  much  more  consistent  and  reacted 
more  like  the  4NQ  solution.  All  y-displacements  for  the 
8NQ  structural  mesh  (coarse)  were  downward  after  680  time 
steps  (t  *  .68E-3). 

The  fourth  variation  undertaken  to  further  improve  the  results 

between  the  4NQ  and  8NQ  solutions  was  new  discretizations  of  the 

structural  meshes  using  more  nodes  and  elements.  These  new 

discretizations  were  done  in  order  to  determine  if  the  element  sizes  in 
the  two  initial  structural  meshes  caused  the  discrepancies  in  the 
results  shown  previously.  The  first  new  discretization  involved  an 

increase  in  the  number  of  nodes  and  elements  in  the  8NQ  structural  mesh 
(the  finer  mesh  in  Figure  3.22b).  An  analysis  using  this  finer  mesh  and 
the  1.0  psi  peak  load  was  executed.  The  outcome  from  this  analysis  was 
much  more  consistent  (displacements  [stresses]  within  the  structure  were 
approximately  equal  at  similar  levels  in  the  mesh)  compared  to  prior 
analyses.  Qualitatively,  this  solution  exhibited  the  same  trends  when 
compared  to  the  4NQ  solution  for  the  1.0  psi  load  performed  previously. 
When  these  two  solutions  were  compared  numerically  at  the  same  time 
step,  the  values  for  the  y-displacements  for  the  8NQ  solution  were  an 
average  of  20"  less  than  those  for  the  4NQ  solution  (stress  values  were 
correspondingly  less  also).  Table  3.6  shows  values  for  the 
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TaMe  3.6  Comparison  of  Y-Displacements  Between  the  8NQ  Solution 
(Finer  Mesh)  and  the  4NQ  Solution  for  a  1.0  psi  Peak 
Load. 


Original 

Node  Location 

x-coordinate  y-coordinate 
(in)  (in) 

8NQ 

Y-Displacement 
t=.8xl0"3  sec 
(in  x  10‘4) 

4NQ 

Y-Di splacement 
t=.8xl0"3  sec 
(in  x  10‘4) 

Percent 
Difference 
Compared  to 
4NQ  at 

t=.8xl0‘3  sec 

0.0 

72.0 

-.1927 

f 

r\> 

►— » 
r\> 

8.8% 

3.0 

72.0 

-.1851 

... 

6.0 

72.0 

-.1749 

-.2009 

12.9% 

9.0 

72.0 

-.1633 

... 

12.0 

72.0 

-.1678 

-.1972 

14.9% 

0.0 

66.0 

-.1774 

-.2073 

14.4% 

6.0 

66.0 

-.1699 

-.1961 

13.4% 

12.0 

66.0 

-.1587 

-.1912 

17.0% 

0.0 

60.0 

-.1618 

-.1996 

18.9% 

3.0 

60.0 

-.1626 

... 

... 

6.0 

60.0 

-.1593 

-.1878 

15.2% 

9.0 

60.0 

-.1550 

... 

•  •  » 

12.0 

60.0 

-.1463 

-.1831 

20.1% 

0.0 

54.0 

-.1498 

-.1949 

23.1% 

6.0 

54.0 

-.1421 

-.1785 

20.4% 

12.0 

54.0 

-.1304 

-.1710 

23.7% 

0.0 

48.0 

-.1343 

-.1898 

29.2% 

3.0 

48.0 

-.1264 

— - 

— 

6.0 

48.0 

-.1218 

-.1711 

28.8% 

9.0 

48.0 

-.1183 

— 

... 

12.0 

48.0 

-.1161 

-.1612 

28.0% 

0.0 

42.0 

-.0944 

-.1788 

47.2% 

6.0 

42.0 

-.1095 

-.1639 

33.2% 

12.0 

42.0 

-.1162 

-.1563 

25.7% 

y-displacements  obtained  from  the  two  solutions  at  time  t  *  .8  x  10”J 
sec  (the  last  time  step).  A  comparison  of  the  two  sets  of  values  from 
the  two  solutions  resulted  in  a  percent  difference  ranging  between  9% 
and  47%  (average  was  approximately  20%).  The  reason  for  this 
significant  difference  between  the  two  solutions  for  such  a  small 
loading  was  not  readily  apparent  to  the  author  at  that  point  in  time. 


However,  the  reason  was  determined  later  in  this  investigation.  Further 
refinements  were  made  to  both  the  4NQ  anc  8NQ  structural  meshes  (the 
finer  4NQ  mesh  and  the  finest  8NQ  mesh  in  Figure  3.22).  These  meshes 
were  then  used  in  the  analyses.  The  results  from  these  new  analyses 
still  produced  a  20%  difference  between  the  two  solutions.  The  two 
finest  meshes  were  then  subjected  to  the  initially  applied  blast  load 
(50  ksi  peak  pressure)  to  determine  if  a  20%  difference  between  the  two 
solutions  would  still  occur.  The  comparison  of  the  results  from  these 
two  executions  yielded  a  15%  difference  between  the  same  values  of 
y-displacement.  It  should  be  noted  that  the  stress  values  could  not  be 
compared  at  similar  time  steps  due  to  the  differences  in  the  wave  speeds 
in  the  two  solutions.  One  further  observation  involving  only  the  8NQ 
solution  was  that  the  results  were  much  more  consistent  compared  to 
previous  analyses  using  the  blast  loading. 

The  fifth  variation  used  In  the  analyses  of  the  8NQ  solution  was  an 
increase  in  the  order  of  integration  from  2.0  to  3.0.  The  integration 
order  for  the  8NQ  was  increased  in  an  effort  to  further  reduce  the 
differences  in  values  between  the  4NQ  and  8NQ  element  solutions.  The 
results  occurring  from  this  change  showed  no  relative  improvement  in  the 
comparison  between  the  two  solutions. 

Three  conclusions  and  one  related  recommendation  were  made  based  on 
the  results  of  the  analyses  for  the  soil-structure  interaction  problem 
using  the  8NQ  element.  The  first  conclusion  was  that  the  8NQ  solutions 
(those  obtained  using  the  two  refined  meshes  in  Figure  3.22b)  produced 
consistent  results  for  the  different  analyses  of  a  soil -structure 
interaction  problem  involving  an  applied  blast  load  and  analyzed 
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axi symmetrical ly  using  the  AFWL  "engineering"  model.  The  results  from 
these  solutions  have  not,  as  yet,  been  shown  to  be  reliable.  In  fact, 
these  solutions  were  significantly  different  (15-20*)  when  compared  to 
those  obtained  with  the  use  of  the  4NQ  element.  These  differences  must 
be  interpreted  carefully  however,  because  the  true  reliability  of  the 
8NQ  solution  has  not  been  verified  since  no  analytical  solutions  or  test 
data  were  available  for  comparison.  The  second  conclusion  was  that  the 
energy  formulation  was  in  error.  The  third  conclusion  was  that  the  8NQ 
element  solutions  were  definitely  obtained  less  efficiently  than 
analyses  obtained  using  the  4NQ  element.  The  8NQ  element  solutions 
required  on  the  average  75%  more  execution  time  than  those  solutions 
using  the  4NQ  elements.  The  reasons  for  this  CPU  time  difference  were 
explained  in  the  discussion  for  the  cantilever  beam.  The  one 
recommendation  was  that  the  8NQ  element  should  be  used  by  personnel  at 
the  AFWL  in  an  analysis  of  a  soil -structure  interaction  problem  for 
which  test  data  are  available.  The  comparison  between  the  test  data  and 
this  8NQ  element  solution  would  give  a  measure  of  the  reliability  of  the 
8NQ  solution.  However,  if  the  8NQ  is  still  found  to  be  unreliable  it 
may  be  due  to  one  of  the  following:  (1)  an  error  in  the  8NQ 
axi symmetric  formulation,  (2)  an  incorrect  interaction  between  the 
slideline  and  8NQ  formulations,  or  (3)  an  incorrect  definition  of 
plastic  flow  in  the  8NQ  formulation.  It  should  be  noted  that  the 
version  of  SAMS0N2  used  in  the  analyses  in  this  study  does  not  include 
any  updates  in  the  slideline  routines  which  have  been  formulated  by  the 
AFWL  in  the  last  couple  years. 
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3.4  Additional  Analyses 

Two  more  problems  were  analyzed  using  the  8NQ  higher-order  element 
(a  4NQ  analysis  was  also  performed  for  each)  for  which  analytical 
solutions  exist.  The  first  problem  was  a  static  analysis  (dynamic 
relaxation)  of  a  fixed-end  beam.  This  beam  had  the  same  dimensions  and 
discretization  as  the  cantilever  beam  discussed  earlier  in  this 
investigation  and  was  subjected  to  a  concentrated  load  at  the  center  of 
the  beam.  The  beam  was  analyzed  using  a  biaxial  elastic-perfectly 
plastic  plane  stress  (Poisson's  ratio  *  0.0)  material  law.  This  problem 
was  chosen  for  analysis  in  order  to  further  determine  the  reliability  of 
a  8NQ  solution  after  plastic  flow  had  been  initiated.  The  sec. id 
problem  was  an  axisynmetric  analysis  involving  one-dimensional  wave 
propagation.  The  orientation  of  the  geometric  configuration  for  the 
previous  one-dimensional  wave-propagation  problem  was  rotated  90°  for 
use  in  the  current  axisymmetric  analysis.  The  same  element 
discretizations  (nodes  renumbered)  and  material  law  were  used  for  this 
axisymmetric  analysis.  The  same  displacement  function  versus  time  was 
applied  to  the  lower  nodes  for  this  axisynmetric  configuration.  The 
boundary  conditions  were  changed  to  fixed  in  the  x-direction  and  free  in 
the  y-direction.  This  problem  was  selected  in  order  to  further  test  the 
reliability  of  a  dynamic  axisymmetric  8NQ  solution.  The  input  data  used 
in  the  analyses  of  these  two  problems  are  in  Appendix  A. 

The  results  from  the  8NQ  element  analysis  of  the  fixed-end  beam 
compared  very  well,  both  elastically  and  plastically,  U  the  analytical 
solutions.  The  maximum  y-displacement  and  normal  stress  values  obtained 
from  the  elastic  8NQ  solution  differed  by  0.5%  and  1.5%,  respectively. 


when  compared  to  values  obtained  from  calculations  involving 
Equations  3.7  and  3.8. 


PL 0  .  1.2  PL 

''max  =  TS2FT  + 


(3.7) 


where. 


peak  value  of  the  applied  concentrated  load  (lbs), 

length  of  the  beam  (ft), 

modulus  of  elasticity  for  the  beam  (psf), 

4 

moment  of  inertia  about  the  x-axis  (ft  ), 
shearing  modulus  (psf),  and 

2 

cross-sectional  area  of  the  beam  (ft  ). 


V  .  /4Px  PL, 
p —  (~g — 


(3.8) 


where. 


*  moment  at  any  distance  x  *  *■  from  the  left  side  of  the 
beam  (ft-lb), 

=  distance  to  the  right  of  the  left  side  of  the  beam  (ft), 
and 

=  distance  above  (+)  or  below  (-)  the  neutral  axis  of  the 
beam  (ft). 


Other  values  for  other  locations  within  the  beam  also  compared  well. 

The  elastic-perfectly  plastic  8NQ  solution  was  not  compared 
directly  to  analytical  solution  values,  but  instead  was  compared 
Qualitatively  to  corresponding  theoretical  concepts  and  to  the  4NQ 
solution  through  observation  of  the  stress  patterns  in  the  8NQ  solutions 
after  the  yield  point  value  had  been  reached.  The  4NQ  and  8NQ  solutions 
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both  exhibited  local  yielding  for  elements  in  their  respective  meshes 
after  the  yield  point  value  had  been  exceeded.  The  stress  (o  ) 
distributions  for  particular  cross-sections  of  the  beam  were  plotted 
both  before  and  after  yielding  and  the  results  were  compatible  with 
mechanics  of  materials  principles  (7,  11).  The  distributions  for  the 
elastic  cross-sections  varied  linearly,  while  those  for  the 
elastic-perfectly  plastic  solutions  were  constant  for  the  top  and  bottom 
portions  of  the  beam  with  connecting  linear  variations.  Another 
observation  made  from  the  results  was  that  redistribution  in  elemental 
stresses  occurred  as  a  result  of  local  yielding.  The  elements  nearest 
to  the  fixed-ends  and  nearest  to  the  center  of  the  beam  had  very  small 
increases  in  stress  values  after  the  initial  elastic  load  was  increased, 
while  the  other  elements  had  much  larger  stress  increases.  This 
difference  between  the  relative  increases  in  stress  values  between  the 
different  elements  was  compatible  with  plasticity  theory.  Two  other 
observations  were  made:  (1)  the  y-displacements  for  the  nodes  at  the 
center  of  the  beam  for  both  the  4NQ  and  8NQ  solutions  were  much  larger 
after  yielding  had  occurred  as  opposed  to  the  elastic  displacements,  and 
(2)  the  output  values  (displacements  and  stresses)  were  larger  in  the 
8NQ  solution  compared  to  those  of  the  4NQ  due  to  the  fact  that  the  8NQ 
solution  predicts  higher  elastic  stress  values  because  of  the  difference 
in  the  order  of  integration  between  the  4NQ  and  8NQ  elements. 

The  results  for  the  axisymmetric  analysis  of  a  wave  propagation 
problem  were  compared  numerically  and  qualitatively  to  the  previous 
solutions  (Section  3.1).  Output  values  obtained  from  the  current  4NQ 
solution  (displacements,  velocities,  and  stresses  in  the  y-direction) 
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were  identical  to  the  corresponding  values  from  the  previous  4NQ 
analysis  of  the  one-dimensional  wave  propagation  problem  (displacements, 
velocities,  and  stresses  in  the  y-direction) .  However,  similar  values 
obtained  from  the  current  axisymmetric  8NQ  solution  did  not  compare  well 
with  those  from  either  the  current  axisymmetric  4NQ  solution  or  the 
previous  one-dimensional  solution.  Figures  3.23,  3.24,  and  3.25  show 
representative  examples  of  the  deviations  between  the  4NQ  and  8NQ 
axisymmetric  solutions.  Similar  deviations  occurred  between  the  two 
solutions  for  other  nodal  y-displacement  and  y-velocity  values  as  well 
as  for  other  elemental  stress  (o  )  values.  In  comparing 
Figures  3.23-3.25  to  similar  figures  in  Section  3.1,  it  was  apparent 
that  the  8NQ  solution  lagged  the  4NQ  solution  in  the  axi symmetrical 
analyses,  particularly  in  the  displacement  and  stress  plots.  The 
results  from  comparing  these  two  axisymmetric  solutions  numerically  at 
the  same  time  step  were  that  the  8NQ  values  were  around  15%  less  than 
those  from  the  4NQ  solution. 

Two  conclusions  were  made  based  on  these  results.  The  first 
conclusion  was  that  the  8NQ  performed  properly  when  used  in  a  static 
elastic-perfectly  plastic  plane  stress  solution.  Strain  hardening  was 
not  attempted  since  the  elastic-perfectly  plastic  solution  results  were 
so  good.  The  second  conclusion  was  that  an  error  appears  to  be  present 
in  the  8NQ  axisymmetric  formulation.  This  conclusion  was  made  on  the 
basis  of  results  obtained  from  the  current  wave  propagation  analyses  and 
the  axisymmetric  analyses  oV  the  soil-structure  interaction  problem. 
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3.5  Summary  of  Results  and  Conclusions 

Five  different  problems  of  various  complexities  were  analyzed  using 
the  current  8NQ  higher-order  element  formulation  in  the  SAMS0N2  code. 
The  5NT  and  6NT  higher-order  element  formulations  were  also  tested  with 
a  cantilever  beam  solution. 

The  results  from  the  analyses  using  the  5NT,  6NT  and  8NQ 
higher-order  elements  were  quite  good  overall.  The  8NQ  element  results 
compared  exceptionally  well  with  analytical  solution  values  for  the 
dynamic  elastic  analysis  of  a  one-dimensional  wave  propagation  problem, 
the  static  and  dynamic  elastic  analyses  of  a  cantilever  beam  problem, 
and  the  static  elastic-perfectly  plastic  analysis  of  a  fixed-end  beam 
problem.  Discrepancies  in  the  8NQ  results  were  found  in  the  dynamic 
axisymmetric  elastic-plastic  analysis  of  the  complex  soil-structure 
interaction  problem  and  the  dynamic  axisymmetric  elastic  analysis  of  a 
one-dimensional  wave  propagation  problem.  The  6NT  element  results 
compared  very  well  with  values  obtained  from  the  analytical  solutions 
for  both  the  static  and  dynamic  elastic  analyses  of  the  cantilever  beam 
problem.  The  5NT  element  results  compared  poorly,  in  general,  to 
analytical  solution  values  for  the  cantilever  beam  analyses.  However, 
some  5NT  values  did  compare  well.  In  particular,  the  maximum  value  for 
the  y-displacement  obtained  for  the  5NT  from  the  dynamic  analysis  of  the 
cantilever  beam  compared  the  best  (compared  to  the  values  from  the  6NT 
and  8NQ  analyses)  to  the  value  calculated  for  an  analytical  solution 
which  neglected  shear  reformation,  axial,  and  rotary  inertia  effects. 
Further  results  showed  that  the  solutions  involving  the  higher-order 


elements  were  less  efficient  compared  to  corresponding  4NQ  element 
solutions. 

Some  difficulties  were  encountered  prior  to  and  during  the 
execution  of  the  higher-order  element  analyses.  The  following  is  a  list 
of  some  of  these  difficulties: 

(1)  The  mesh  generation  routine  in  SAMS0N2  was  not  formulated 
so  that  the  meshes  used  in  the  higher-order  element 
solutions  could  be  generated  easily.  For  instance,  the 
8NQ  elements  can  only  be  generated  perpendicular  to  the 
direction  in  which  the  nodes  were  generated  due  to  the 
fact  that  these  elements  have  different  node  increments 
between  corner  nodes  and  midside  nodes.  Another  example, 
was  that  the  5NT  element  mesh  used  in  the  cantilever  beam 
analyses  could  not  be  generated  for  reasons  similar  to 
those  stated  for  the  8NQ  element. 

(2)  The  placement  of  the  midside  nodes  was  critical  in  order 
to  obtain  a  solution  using  the  5NT  elements.  It  was 
mentioned  previously  in  Section  3.2  that  an  attempt  was 
made  to  improve  a  5NT  solution  by  changing  the  locations 
of  the  midside  nodes.  The  elemental  nodes  were  input 
according  to  Reference  3  for  these  rearranged  nodes,  but 
negative  areas  were  calculated  for  the  elements.  This 
calculation  occurred  during  the  nodal  mass  allocation, 
and  therefore  was  not  affected  by  *ue  time  step  size.  A 
solution  was  never  obtained  for  the  mesh  containing  the 


reordered  nodes. 
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(3)  In  the  execution  phase  of  each  analysis,  it  was  observed 
that  the  higher-order  element  analyses  required  a  smaller 
time  step  size  as  compared  to  the  4NC  element  analyses. 

The  time  step  value  used  for  the  cantilever  beam  analyses 
was  set  equal  to  the  maximum  value  required  to  obtain  a 
stable  solution  using  the  8NQ  as  determined  by  a  trial 
and  error  procedure.  The  stable  time  step  for  the  4NQ 
solution  was  larger  even  though  the  4NQ  mesh  used  in  the 
analyses  was  significantly  finer. 

A  summary  of  the  conclusions  discussed  earlier  in  this  chapter  is 
presented  here.  The  higher-order  element  (5NT,  6NT,  and  8NQ) 
formulation  for  the  calculation  of  strain  appears  to  be  correct  for  both 
static  and  dynamic  plane  analyses.  The  8NQ  element  formulation  for  the 
internal  force  calculations  (stress  calculations  are  part  of  the 
internal  force  calculations)  was  shown  to  be  functioning  properly  for 
static  analyses  using  the  biaxial  elastic-plastic  plane  stress  material 
law  and  for  dynamic  analyses  using  the  biaxial  elastic  plane  strain  and 
plane  stress  material  laws.  Internal  force  calculations  were  also 
executed  correctly  for  the  6NT  element  in  static  and  dynamic  analyses 
using  the  biaxial  elastic  plane  stress  material  law.  These  calculations 
were  also  correct  for  the  5NT  element  assuming  that  the  5NT  element 
cantilever  beam  solution  could  have  been  improved  by  increasing  the 
number  of  elements  in  the  mesh.  The  8NC  axisymmetric  formulation  was 
found  to  be  incorrect  for  dynamic  analyses  using  the  AFWL  "engineering" 
model  and  the  elastic  plane  strain  material  laws.  Three  potential 


additional  errors  were  discovered  in  the  SAMS0N2  finite  element 
formulation  and  are  as  follows: 

(1)  The  energy  error  formulation  is  incorrect  in  its  present 
state  as  was  shown  in  the  analysis  of  the  soil 
structure-interaction  problem. 

(2)  The  slideline  calculations  for  the  8NQ  analysis  appeared 
to  be  quite  different  when  compared  to  those  obtained  for 
the  4NQ  analyses. 

(3)  The  AFWI.  "engineering"  model  might  also  be  in  error  for 
the  8NQ  formulation  when  used  for  elastic-plastic 
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CHAPTER  4 


Verification  of  the  SAMSON2  Finite  Element  Formulation  for  the 


Eiaht-Node  Quadri 1 ateral  Isoparametric  Continuum  Element 


The  contents  of  this  chapter  dealing  with  the  verification  of  the 
finite  element  formulation  for  the  8NQ  isoparametric  continuum  element 
are  separated  into  three  main  sections.  The  first  section  provides  the 
necessary  background  information  needed  to  better  understand  the 
capabilities  of  the  SAMS0N2  code.  This  information  was  also  important 
in  the  code  verification.  The  second  section  contains  a  synopsis  of  the 
equations  involved  in  the  8NQ  finite  element  formulation  currently  in 
the  SAMS0N2  code.  Each  equation  presented  was  compared  with  information 
existing  in  finite  element  texts,  journal  articles,  and  other  sources  in 
order  to  verify  the  formulation  (1,  5,  7,  8,  9,  11-18).  Any 
disagreements  obtained  from  these  comparisons  are  pointed  out  and  the 
conflicting  equations  from  the  reference  materials  are  presented. 
Solution  flow  charts  and  some  discussion  of  the  equations  are  also 
introduced  in  this  section.  The  third  section  in  this  chapter  includes 
a  summary  of  all  of  the  errors  found  during  the  verification  of  the  8NQ 
finite  element  formulation. 


4.1  Background  Information 

The  SAMS0N2  code  is  a  dynamic  nonlinear  two-dimensional 
structu, e-media  interaction  computer  code  (i)  which  uses  an  explicit 
central  difference  finite  element  solution  scheme.  Explicit  integration 
is  only  conditional ly  stable.  Therefore,  an  accurate  solution  can  only 


be  obtained  by  selecting  a  small  time  step  At.  An  energy  balance  is 
performed  at  the  end  of  each  time  step  to  check  if  the  solution  is 
stable  (energy  error  less  than  1.0%).  A  diagonal  lumped  mass  matrix  is 
used  in  accordance  with  the  explicit  integration  scheme  so  that  no 
solution  of  any  simultaneous  equations  is  required  in  advancing  a  time 
step.  A  stiffness  matrix  is  never  computed  in  an  explicit  time 
integration  scheme  due  to  the  fact  that  solutions  of  simultaneous 
equations  are  never  performed. 

The  SAMS0N2  code  was  designed  specifically  for  analyses  of  large 
displacements,  large  strain  problems  involving  nonlinear  material 
behavior  and  structure-media  interface  (SMI)  problems.  The  stress  and 
strain  calculations  are  performed  in  the  analyses  of  these  nonlinear 
type  problems  with  the  use  of  the  Cauchy  stress  and  velocity  strain 
(rate  of  deformation)  tensors.  The  evaluation  of  stresses  is  performed 
in  a  corotational  coordinate  system  which  rotates  coincident  with  the 
rotation  of  a  quadrature  point  in  an  element.  Therefore,  a  Jaumann  type 
correction  is  not  needed  to  maintain  objectivity  in  the  solution.  The  x 
and  y  nodal  coordinates  used  In  the  analyses  are  the  spatial  Eulerian 
coordinates  which  are  consistent  with  the  velocity  strain  tensor. 
Gaussian  quadrature  is  used  in  the  numerical  integration  of  the  finite 
element  equations  for  the  stress  and  strain  tensors  for  the  8NQ 
isoparametric  continuum  element. 

Many  other  features  are  present  in  SAMS0N2  that  are  linked  to  the 
general  finite  element  formulation  in  the  code.  The  following  is  a  list 
of  some  of  these  additional  features: 


(1)  Multiple  time  step  integration  is  available  for  problems 
containing  more  than  one  mesh  size  in  order  to  eliminate 


unnecessary  integration  of  the  coarser  meshes. 


(2)  Dynamic  relaxation  techniques  are  employed  in  the  SAMS0N2 


code  in  order  to  reduce  the  dynamic  analyses,  with  the 


use  of  sufficient  damping,  to  a  static  equilibrium 


solution. 


(3)  Mass  and  stiffness  (artificial  viscosity)  proportional 


damping  are  used  in  the  SAMS0N2  code. 


(4)  A  single  large  storage  array  (Q  array)  is  used  to  store 


most  of  the  major  variable  arrays  used  in  the  execution 


of  SAMS0N2. 


(5)  SAMS0N2  contains  many  material  models  which  range  from 


the  relatively  simple  elastic-plastic  uniaxial  stress 


material  law  to  the  complicated  viscoplastic  material 


(6)  Slideline  interface  routines  are  used  in  the  SAMS0N2  code 


to  model  relative  sliding  motion  or  material  separation 


between  two  different  types  of  materials  such  as  occurs 


in  SMI  problems.  These  interfaces  can  also  be  rigid  in 


order  to  connect  two  different  type  elements. 


4.2  8NQ  Isoparametric  Continuum  Element  Formulation  and  Verification 


This  section  contains  the  documentation  of  the  research  which  was 


performed  during  the  verification  of  the  current  SAMS0N2  finite  element 


formulation  for  the  8NQ  isoparametric  continuum  element.  The  following 


a 


£ 


is  a  list  of  the  particular  items  which  were  investigated  in  the 
formulation: 


(1)  the  equations  involved  in  the  determination  of  the  total 
diagonal  lumped  mass  matrix,  [M], 

(2)  the  equations  involved  in  the  determination  of  the 
velocity  strains,  (dl,  Cauchy  stresses , (a} ,  and  internal 
forces,  {Fint>, 

(3)  the  calculations  involved  in  the  determination  of  the 
nodal  forces  due  to  the  externally  applied  loads,  (F  }, 

(4)  the  mass  and  stiffness  proportional  damping  (C^M]  + 
C2[K])  calculations, 

(5)  the  equations  used  in  calculating  the  internal  strain 
energy  (U )  and  the  external  work  (W)  terms,  and 

(6)  the  solution  of  the  equation  of  motion  for  the  nodal 

••  • 

accelerations  tub  velocities  (ul,  and  displacements  {u > . 

These  items  were  investigated  in  detail  and  are  discussed  In  the 
following  subsections. 

4.2.1  Element  Configuration  and  Shape  Functions 
and  the  Derivatives  for  the  8NQ 

Some  Important  information  relevant  to  many  phases  of  the 
formulation  development  for  the  8NQ  isoparametric  continuum  element  is 
presented  here.  The  configuration  of  the  8NQ  element  for  which  the 
formulation  was  developed  is  shown  in  Figure  4.1.  The  local  (S,n) 
coordinate  system  is  shown  at  the  center  of  the  element.  This 
coordinate  system  is  defined  such  that  the  values  for  i  and  r  vary 


between  1.0  (value  on  the  right  and  upper  element  boundaries)  and  -1.0 
(value  on  the  left  and  lower  element  boundaries). 
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Figure  4.1  8NQ  Isoparametric  Continuum  Element. 

The  shape  (interpolation)  functions  (Nj)  for  the  8NQ  element  are  shown 
in  Table  4.1  in  terms  of  the  local  coordinate  system.  These  shape 
functions  are  used  in  both  the  interpolation  of  the  element  coordinates 
(Equation  4.1)  and  the  element  displacements  (Equation  4.2)  for  an 
isoparametric  finite  element  formulation. 


x,y  =  global  (x,y)  coordinates  for  any  point  in  any  element, 

NT  *  shape  (interpolation)  function  corresponding  to  node  I, 
and 

x j ,y j  =  global  (x,y)  coordinates  for  node  I. 


8 

u  =  I  N.u  . 
x  T  ,  I  xl 
±  •  X 


where. 


Uxl*  Uyl 


uy  =  Niuyi 


(4.2) 


x  and  y  displacements  for  any  point  in  an  element. 


u  ,u  * 

x  y 

N,  =  shape  (interpolation)  function  corresponding  to  node 

1  I,  and 


x  and  y  displacements  for  node  I. 


The  fundamental  property  associated  with  these  shape  functions  is  that 
the  value  for  is  unity  at  node  I  and  is  zero  at  all  other  nodes.  The 
derivatives  of  these  shape  functions  with  respect  to  the  local  element 
coordinates  (3Nj/H  and  3Nj/3n)  are  also  shown  in  Table  4.1.  These 
derivatives  are  provided  because  they  are  more  important  than  the  shape 
functions  in  the  actual  finite  element  formulation.  As  a  final  note  to 
this  section,  the  shape  functions  and  their  derivatives  from  Table  4.1 
exist  in  the  SAMS0N2  code  correctly  and  were  verified  through 
comparisons  with  References  1  and  12. 


4.2.2  Determination  of  the  Total  Diagonal  Lumped  Mass  Matrix 


This  section  gives  the  equations  used  i..  'jAMS0N2  to  determine  the 
total  diagonal  lumped  mass  matrix,  [Ml,  for  the  8NQ  element.  Table  4.2 
provides  a  flow  chart  of  the  main  subroutines  involved  in  the 


calculation  of  [Ml  and  the  operations  performed  in  each  subroutine.  The 
mass  matrix  [M]  is  developed  only  once  and  is  assembled  prior  to  any 
other  calculations  in  the  SAMS0N2  solution  scheme.  Therefore,  the 
elemental  masses  are  determined  from  the  initial  problem  configuration. 


Table  4.2  SAMS0N2  Flow  Chart  for  the  Determination  of  the  Total 
Diagonal  Lumped  Mass  Matrix  [Ml  for  an  8NQ  Isoparametric 
Plane  or  Axisymmetric  Continuum  Element. 


Subroutine 


Operation 


1.  ASSBLE  (called  by  SMAIN) 
DO  I  =  1,  #  OF  ELEMENTS 


2.  VASME 


3.  V8ASME 


4.  GAUSS1 


CONTINUE 


Assembles  the  elemental  mass  matrices 
into  the  total  diagonal  lumped  mass 
matrix  [M]. 

This  routine,  which  is  called  by  ASSBLE, 
calls  the  appropriate  VnASME  routine 
which  is  used  to  compute  the  mass  matrix 
for  an  element  which  contains  n  nodes. 

This  routine  is  used  to  compute  the 
elemental  mass  matrix  for  an  8NQ 
isoparametric  plane  or  axisymnetric 
continuum  element  for  a  particular  order 
of  integration  (iorder). 

This  routine  provides  the  Gauss-Legendre 
abscissae  and  weight  coefficient  values 
based  on  Iorder  which  are  used  in  V8ASME 
in  order  to  determine  the  elemental  mass 
matrices. 

1=1+1. 


The  two  subroutines  of  primary  importance  in  this  investigation 
were  V8ASME  and  GAUSS1.  The  subroutine  V8ASME  contains  all  of  the 
formulation  used  in  the  calculation  of  [M],  GAUSS1  provides  the 
necessary  Gauss-Legendre  abscissae  and  weight  coefficients  used  in  the 


numerical  integration  of  [M].  These  coefficients  are  shown  in 
Appendix  B  with  the  column  labeled  ia  containing  the  abscissae 
coefficients  and  the  column  labeled  h  containing  the  weight 
coefficients.  These  coefficients  exist  in  the  subroutine  GAUSS1  to  15 
decimal  places.  However,  some  of  these  coefficients  were  incorrectly 
typed  into  the  SAMS0N2  code.  Table  4.3  shows  a  list  of  the  coefficients 
that  were  found  to  be  incorrect. 

The  first  step  in  the  process  of  determining  [M]  was  to  calculate 
the  mass  contained  in  each  element  (performed  in  V8ASME). 
Equations  4.3,  4.4,  and  4.5  were  used  to  calculate  the  area  of  an 
element. 


i order  i order 

A  =  i  E  h  h.  | J | 
i«l  j=l  J 


(4.3) 

(4.4) 


3x  _ 

8  3N. 

V  ^  v 

Li  , 

8 

3N  j 

r 

H 

i=i  3<  1 

3n 

L 

1*1 

3n 

3  V  . 

8  3N. 

E  — -  V, 

ll  - 

8 

3  Nj 

z 

H 

i.i95  1 

3n 

1*1 

3n 

where, 

A 

i order 


(4.5) 


area  of  an  element, 

order  of  integration  selected  for  the  problem 
analyses, 

Gauss-Legendre  weight  coefficients, 

determinant  of  the  Jacobian  matrix  which  contains 
the  derivatives  of  the  global  (x,y)  coordinates  with 
respect  to  the  local  (s '  coordinates. 
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Table  4.3  List  of  Gauss-Legendre  Abscissae  and  Weight  Coefficients 
Incorrectly  Typed  into  the  Subroutine  GAUSS1. 


Order  of 
Integration 
n 

SAMS0N2 

Coefficient 

Correct 

Coefficient 

2 

-0.57735 

02691 

8926 

-0.57735 

02691 

89626 

4 

+0.33998 

81043 

584856 

+0.33998 

10435 

84856 

8 

+0.96202 

89856 

497536 

+0.96028 

98564 

97536 

9 

+0.83603 

11073 

2663 

+0.83603 

11073 

26636 

9 

+0.96816 

02395 

07656 

+0.96816 

02395 

07626 

10 

-0.67940 

65682 

99024 

-0.67940 

95682 

99024 

a)  Abscissae 

i  Coefficients 

Order  of 

SAMS0N2 

Correct 

Integration 

n 

Coefficient 

Coefficient 

3 

0.55555 

55555 

556 

0.55555 

55555 

55556 

3 

0.88888 

88888 

886889 

0.88888 

88888 

88889 

3 

0.55555 

55555 

556 

0.55555 

55555 

55556 

8 

0.22238 

10345 

3374 

0.22238 

10344 

53374 

10 

0.26926 

68193 

09996 

0.26926 

67193 

09996 

10 

0.24908 

63625 

15982 

0.21908 

63625 

15982 

b)  Weight  Coefficients 


v.v 


derivatives  of  the  global  (x,y)  coordinates 
with  respect  to  the  local  (C,n)  coordinates, 


11  11  11  11 
3C’  3 r  *  3C*  3n 


3Nj  SNj 


shape  function  derivatives  *rom  Table  4.1,  and 


xT,y.  =  x  and  y  coordinates  for  any  node  I  in  an  8NQ 
element. 


Equations  4.4  and  4.5  were  evaluated  using  the  Gauss-Legendre  abscissae 
coefficients  obtained  from  GAUSS1  for  a  particular  order  of  integration 
(iorder).  These  coefficients  were  substituted  into  these  equations  for 
C  and  n  for  all  possible  combinations  of  the  coefficients.  Equation  4.3 
requires  the  use  of  both  the  Gauss-Legendre  abscissae  and  weight 
coefficients  in  order  to  numerically  integrate  the  area  of  an  element. 
These  three  equations  used  in  the  SAMS0N2  formulation  were  verified 
through  comparisons  with  References  14  and  18.  The  volume  and  the 
corresponding  mass  of  an  element  were  computed  upon  completion  of  the 
area  calculation.  The  volume  calculation  for  a  plane  element  involved 
multiplying  the  area  by  the  thickness.  The  imss  of  an  element  was  then 
calculated  by  multiplying  the  volume  times  the  mass  density  as  shown  in 
Equation  4.6. 

iorder  iorder 

M  »  ot  l  l  h.h.  j  J  |  (4.6) 

i*l  j-1  1  J 

where, 

M  =  mass  for  a  plane  element, 

p  =  mass  density  of  the  element,  and 

t  =  thickness  of  the  element. 


The  volume  calculation  for  an  axisymnetric  element  involved  multiplying 
the  area  by  2r.  times  the  distance  from  tne  axis  of  symmetry  to  the 


geometric  center  of  an  element  (r).  However,  the  SAMS0N2  formulation 
for  determining  r  was  found  to  be  incorrect.  The  correct  formulation 
for  the  determination  of  r  is  shown  in  Equation  4.7. 


r  = 

where , 
r 


X1  +  x2  +  X3  +  x4 

- o - 


distance  from  the  axis  of  symmetry  to  the  geometric 
center  of  an  8NQ  element,  and 


(4.7) 


;.,x-,x-,x.  =  x-coordinates  for  the  four  corner  nodes  of  an  8NQ 
1  element  (Figure  4.1). 


The  SAMS0N2  formulation  used  a  value  of  3.0  instead  of  the  value  4.0  in 
the  denominator.  This  error  was  verified  by  comparisons  to  formulae  in 
References  1  and  14  and  through  comparisons  with  hand  calculated  values. 
In  addition  to  this  error,  the  factor  2*  was  not  used  in  the  calculation 
of  the  volume.  However  the  elimination  of  2*  was  justified  in  that  it 
canceled  out  in  the  final  solution  of  the  equation  of  motion.  A  further 
note  regarding  the  ax i symmetric  volume  calculation  is  that  the 
calculation  using  r  is  only  an  approximation  of  the  actual  value.  The 
mass  was  computed  for  an  axisymmetric  element  using  Equation  4.8.  This 
equation  with  the  use  of  the  correct  r  value  was  verified  as  being 
correct. 

i order  i order 

M=o?E  E  h.h.jjj  (4.81 

i=l  j=l  1  J 


The  second  step  in  the  process  of  determining  [Ml  was  the 
allocation  of  the  calculated  element  mass  to  the  appropriate  nodes.  ">e 
element  mass  was  allocated  to  the  nodes  of  an  8NQ  according  v 
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Figure  4.2.  As  illustrated,  20?  and  5 %  of  the  element  mass  were 
allocated  to  the  midside  nodes  and  corner  nodes,  respectively.  This 
nodal  mass  distribution  used  in  the  SAMS0N2  coce  could  not  be  verified. 
However,  it  was  found  to  be  similar  to  a  tributary  area  approach  where 
the  values  were  18.75?  and  6.25?.  It  was  also  *ound  to  be  similar  to  an 
approach  which  used  the  scaled  diagonal  terms  of  a  consistent  mass 
matrix  where  the  values  were  22.22?  and  2.78?.  In  reviewing  the 
different  references  (5,  9,  15),  it  was  determined  that  no  standard 
lumping  scheme  for  the  8NQ  existed.  Instead,  many  schemes  were  found 
that  provided  good  results  when  used  in  problem  solutions  involving  the 
8NQ  continuum  element. 


.05  .20  .05 
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Figure  4.2  Nodal  Lumped  Masses  for  an  8NQ  Element. 

The  final  step  in  the  process  of  determining  [M]  was  the  storage  of 
these  nodal  masses.  The  mass  allocated  to  each  node  in  an  element  was 
stored  for  each  nodal  degree  of  freedom  in  an  array  called  smass.  The 
nodal  mass  allocations  determined  in  step  twc  were  added  to  smass  for 
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each  element.  The  diagonal  lumped  mass  matrix  [M]  was  completed  when 
the  mass  contribution  of  every  element  was  stored  in  the  smass  array. 
Upon  the  determination  of  [M],  this  mass  matrix  remained  unchanged 
throughout  the  entire  solution  and  was  only  used  to  solve  the  equation 
of  motion  at  the  end  of  each  time  step. 


4.2.3  Determination  of  the  Velocity  Strains 

This  section  displays  the  equations  formulated  in  SAMS0N2  that  are 
used  to  calculate  the  velocity  strain  tensor,  (d).  The  velocity  strain 
tensor  measures  the  current  rate  of  deformation  and  is  used  for  problems 
involving  geometric  nonlinearities  (large  rigid  body  rotations  and 
deformations).  Table  4.4  shows  a  flow  chart  of  the  pertinent 
subroutines  and  their  functions  that  are  used  in  SAMS0N2  in  order  to 
determine  the  internal  nodal  forces  (F^L  The  velocity  strains  (d) 
are  computed  in  subroutine  V8FRCN  in  step  a  as  shown  in  Table  4.4. 

The  general  equations  for  the  velocity  strains  are  similar  to  the 
strain-displacement  relations  for  small -displacement  theory. 
Equation  4.9  shows  the  four  velocity  strain  components  for  an  8NQ 
axisymmetric  continuum  element. 
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Table  4.4  Flow  Chart  tor  the  Determination  of  the  Internal  Nodal 
Forces  and  Internal  Strain  Energy  for  an  8NQ 
Isoparametric  Plane  or  Axisynmetric  Continuum  Element. 


Subroutine 


Operation 


1.  FRCIN  (called  by  SMAIN) 
DO  I  =  1,  #  OF  ELEMENTS 


2.  VFRCIN 


3.  V8FRCN 


Combines  the  internal  nodal  forces  for 
each  element  into  a  total  internal  nodal 
force  vector.  Combines  the  internal 
strain  energies  for  each  element  into  one 
value  for  each  element  type. 

Calls  the  appropriate  element  routine 
(VnFRCN)  which  comoutes  the  internal 
noclal  forces  for  an  element  with  n  nodes. 

Computes  the  internal  nodal  forces  for  an 
8NQ  isoparametric  plane  or  axi symmetric 
continuum  element  using  the  following 
steps: 

(a)  Computes  the  global  components  for  the 
velocity  strains  using  subroutines 
GAUSS 1,  V8BMAT ,  V8SHAP,  and  GMPRD. 

(b)  Computes  the  corotatlonal  strain 
increments  using  subroutines  VSTRAN 
and  GMPRD. 

(c)  Calls  the  appropriate  stress  routines 
through  STRES  which  compute  the  total 
stresses  based  on  the  current  strain 
increments  and  the  viscous  stresses 
due  to  damping  (SDAMP).  Also,  the 
total  strains  are  updated. 

(d)  Computes  the  internal  nodal  forces  and 
the  internal  strain  energy  for  an  8NQ 
element  using  the  stress  values  from 
(c)  and  subroutines  STPRD  and  GTPRD. 


CONTINUE 


1=1+1. 
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where. 


d  ,d  ,d 
x’  y*  z 


Vx’Vy 


VxI’VyI 
3Nj  3Nj 

Tx~'Jy- 


three  normal  velocity  strain  components  in  the  x,  y, 
and  z  (r,  z,  and  e)  directions, 

shearing  velocity  strain  component, 

velocity  components  for  any  point  in  an  element 
(these  velocities  are  never  calculated), 

distance  from  the  axis  of  symmetry  to  the  quadrature 
point  at  which  the  strains  were  computed, 

current  global  deformed  nodal  coordinates  for  an  8NQ 
element, 

velocity  components  for  node  I,  and 


derivatives  of  the  shape  functions  with  respect  to 
the  global  (xf>/  coordinate  system. 


The  velocity  strain  component  In  the  z-direction,  d2,  is  equal  to  zero 
for  analyses  of  plane  continue.  The  derivatives  of  the  shape  functions 
with  respect  to  the  global  coordinate  system  were  determined  using 
Equations  4.10  and  4.11. 
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The  velocity  strain  components  shown  in  Ecuation  4.9  were  evaluated 
using  a  matrix  [B],  This  matrix  incorporates  the  expressions  from 
Equations  4.9,  4.10,  and  4.11  such  that  the  velocity  strain  components 
are  obtained  through  the  multiplication  of  [E]  times  the  nodal  velocity 
components  (VxJ  and  Vyj). 

Table  4.5  shows  the  results  obtained  when  multiplying  the  TB] 
matrix  used  in  SAMS0N2  times  the  nodal  velocities.  As  shown,  the  - 
velocity  strains  are  multiplied  by  a  factor  equal  to  |J|.  The  velocity 
strains  were  not  directly  calculated  because  the  |"B]  matrix  was 
developed  for  use  in  the  determination  of. the  internal  nodal  forces,  so 
that  the  I J I  term  canceled  out  (see  pp.  53-55  in  Reference  1).  The 
determination  of  the  [B]  matrix  was  performed  in  the  V8BMAT  subroutine 
using  values  calculated  for  the  shape  function  derivatives  (3Nj/3x, 
3Nj/3y)  in  subroutine  V8SHAP.  The  [B]  matrix  was  calculated  for  each 
quadrature  point  within  an  element  using  the  Gauss-Legendre  abscissae 
coefficients  obtained  from  the  GAUSS1  routine.  Therefore,  the  velocity 
strains  were  also  computed  for  the  same  quadrature  points.  The  general 
matrix  product  routine  GMPRD  in  SAMS0N2  was  used  to  multiply  the  [Bl 
matrix  by  the  nodal  velocity  matrix  as  shown  in  Table  4.5. 
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The  next  step  performed  in  the  calculation  for  the  velocity  strains 
was  the  transformation  of  these  strains  from  the  global  (x,y)  coordinate 

A  A 

system  to  the  corotational  reference  frame  (x,y).  A  corotational 
reference  frame  exists  at  each  of  the  quadrature  points  and  rotates 
according  to  the  movement  of  these  points.  Therefore,  in  addition  to 
the  velocity  strain  calculations,  a  rotational  velocity  (w)  was  computed 
for  each  quadrature  point  using  Equation  4.12. 


,n+J  _ 
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(4.12) 


Equation  4.12  was  found  to  be  identical  to  the  equation  for  vorticity 
used  in  fluid  mechanics  (18).  The  rotational  velocity  was  multiplied  by 
the  value  for  the  time  step  in  order  to  obtain  the  current  rotation 
increment.  The  total  rotation  for  the  current  time  step  (0n+1)  was  then 
calculated  using  Equation  4.13  and  was  then  stored  in  the  strs  array  in 
the  last  location  assigned  to  an  element  quadrature  point. 


en+1  «  en  +  at  «n+*  (4.13) 

where, 

en+1,en  «  total  rotations  which  have  occurred  at  a  quadrature 
point  for  the  current  and  previous  time  steps, 

at  *  time  step  used  for  the  explicit  integration  solution 

scheme,  and 

wn+*  »  rotational  velocity  for  a  quadrature  point  which 

occurred  during  the  current  time  step. 

The  total  rotation,  on+1,  was  then  used  in  the  transformation  matrix, 
[Rl,  in  order  to  convert  the  velocity  strains  from  the  global  (x,y) 


coordinate  system  to  the  corotational  framework  (x,y).  Equations  4.14 
and  4.15  show  the  strain  transformation  matrix  [R]  and  the  multiplica¬ 
tion  process  which  was  performed  in  order  to  determine  the  corotational 
velocity  strains.  The  [R]  matrix  was  not  multiplied  by  the  velocity 
strain  in  the  z-direction  since  d2  is  not  affected  by  rotation  in  the 
x ,y  plane. 

cos2©  sin2©  cos  ©  sin  0 

[R]  =  sin20  cos2©  -cos  0  sin  ©  (4.14) 

_-2sin  0  cos  0  2s i n  0  cos  ©  cos2©-sin20_ 

(d  |J|)  =  [R]  (d  |J|)  (4.15) 

where, 

(d  | J I >  =  velocity  strains,  multiplied  by  a  factor  |J|,  in  the 

*  A 

(x,y)  coordinate  system, 

[Rl  =  velocity  strain  transformation  matrix,  and 

(d  iJl)  *  velocity  strains,  multiplied  by  a  factor  |J|,  in  the 

global  (x,y)  coordinate  system. 

The  transformation  matrix  was  computed  in  subroutine  VSTRAN  in  the 
SAMS0N2  code.  A  possible  error  was  found  in  this  subroutine  VSTRAN. 
For  values  of  !©!  <  1.0  x  10“5,  the  value  for  each  of  the  sin  0  terms 
was  set  equal  to  zero.  However,  according  to  small  angle  theory, 
sin  0  =  0.  This  error  may  be  trivial  in  the  overall  solution.  The 
multiplication  process  illustrated  in  Equation  4.15  was  executed  in 
subroutine  GMPRD. 

One  last  computation  was  executed  for  the  velocity  strains.  This 
computation  involved  the  multiplication  of  the  corotational  velocity 
strains  by  an  appropriate  factor  so  as  to  obtain  the  incremental  strains 


.'WTrr»".’T%rwv  fv  roww  v  v 


for  the  current  time  step.  This  computation  is  illustrated  in 
Equation  4.16  and  was  the  last  one  performed  prior  to  the  calculation  of 
the  stresses,  {a},  and  the  internal  forces,  }. 


(Ae>  *  {d  IJI } 


(4.16) 


where. 


{Ae}  =  incremental  strains  determined  *cr  the  current  time  step. 

The  formulation  for  the  velocity  strain  calculations  was  found  to 
be  correct  when  compared  to  material  in  Refe-ences  1,  9,  12,  and  16. 
The  only  possible  error  found  in  the  formulation  was  the  use  of  the 
approximation  of  sin  q  *  0.0  for  lej  <  1.0  x  1C~^.  Also,  the  use  of  the 
[B]  matrix  for  the  calculation  of  the  velocity  strains  was  corrected 
later  when  a  factor  of  1/  ! J |  was  multiplied  times  the  velocity  strains. 


4.2.4  Evaluation  of  Stresses  and  Internal  Forces 

The  stresses  calculated  in  SAHS0N2  were  evaluated  using  the  Cauchy 
stress  tensor.  This  tensor  was  computed  in  tne  corotational  coordinate 
system  so  as  to  maintain  objectivity  in  the  overall  solution.  The 
Cauchy  stress  tensor  is  energetically  conjugate  to  the  velocity  strain 
tensor.  Hence,  the  total  internal  work  can  be  calculated  using  these 
tensors.  The  Cauchy  stress  tensor  is  used  fo-  the  incremental  analysis 
of  geometrically  nonlinear  problems. 

The  SAMS0N2  code  contains  eight  differert  constitutive  laws  which 
are  used  to  determine  the  Cauchy  stresses  based  upon  the  incremental 
strains  already  computed.  In  this  study,  only  the  biaxial  material  laws 
(subroutines  STRES3  and  STRES6)  and  the  AFWL  engineering"  material  law 


(subroutine  STRES9)  were  investigated.  The  remainder  of  the  material 
laws  were  not  investigated  because  they  were  considered  less  important 
(uniaxial  laws)  or  had  been  investigated  in  prior  studies.  The 
stress-strain  relations  are  provided  in  the  following  subsections  for 
each  material  law  investigated.  An  example  of  the  derivations  that  were 
required  in  order  to  determine  the  biaxial  elastic-plastic  stress-strain 
matrices  are  also  provided  in  Appendix  C. 

4.2.4. 1  Stress  Evaluations  Using  the  Biaxial  Elastic 
Plane  Strain  Material  Law  (STRES6) 

The  formulation  for  the  STRES6  subroutine  was  developed  for 
isotropic  linear  elastic  materials.  STRES6  can  be  used  for  plane  strain 
or  axisymmetric  analyses.  All  computations  were  performed  in  STRES6 
under  the  assumption  of  totally  elastic  material  behavior.  Therefore, 
no  yield  condition  or  plastic  flow  rule  is  used. 

The  stress-strain  relations  for  the  biaxial  elastic  plane  strain  or 
axisymmetric  material  law  were  expressed  in  terms  of  total  stresses  and 
total  strain  as  opposed  to  Incremental  values  due  to  the  assumption  of 
totally  elastic  behavior.  Hence,  the  current  total  strains  (e}n+*  were 
computed  prior  to  any  stress  calculations  by  adding  the  incremental 
strains  (Ae}n+*  (Equation  4.16)  to  the  total  strains  at  the  end  of  the 
previous  time  step  (e}n  as  shown  in  Equation  4.17. 

(c}n+1  =  {e}n  +  At(d>n+i  (4.17) 


The  total  stresses  were  then  computed  using  the  stress-strain  relations 
shown  in  Equation  4.18. 
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These  equations  were  used  for  both  plane  strain  and  ax i symmetric 
analyses  with  the  strain  in  the  z-direction,  e  ,  being  equal  to  zero  for 
plane  strain  analyses.  The  values  for  e and  were  both  found  to  be 
equal  to  zero  for  plane  strain  analyses  using  v  =  o  (Poisson's  ratio). 
The  current  total  stresses  and  total  strains  for  the  corotational 
coordinate  system  were  stored  in  the  appropriate  locations  within  the 
strs  array  for  each  element  quadrature  point. 

The  calculations  executed  in  STRESS  using  Equations  4.17  and  4.18 
were  all  verified  through  comparisons  to  References  1,  7,  12,  14,  and 
15.  Therefore,  the  biaxial  elastic  plane  strain  or  axisynmetric 
material  law  was  found  to  be  correct. 


4. 2. 4. 2  Stress  Evaluations  Using  the  Biaxial  Elastic-Plastic 


Plane  Stress  Material  Law  (STRES3) 


The  subroutine  STRES3  was  developed  for  analyses  involving 

elastic-perfectly  plastic  or  strain  hardening  materials.  The  solution 

E  EP 

scheme  for  this  subroutine  is  illustrated  in  Table  4.6  with  C  and  C 
referring  to  the  elastic  and  elastic-plastic  stress-strain  matrices 
which  will  be  presented  in  this  section. 


Table  4.6  Solution  Algorithm 
Calcul ations . 


Elastic-Plastic  Stress 


Given:  STRAIN  =  total  strains  for  current  time  step,  U} 

EPS  *  total  strains  for  previous  time  step,  (c}n 

SIG  *  total  stresses  for  previous  time  step,  {o)n 

The  following  procedure  is  used  to  calculate  the  total  stresses  TAU  for  the 
current  time  step,  {o}n+1. 

(a)  Calculate  the  strain  increment  OELEPS,  (ae)n+*: 

DELEPS  «  STRAIN  -  EPS 

(b)  Calculate  the  stress  Increment  DELSIG,  Uo}n+*,  assuming  elastic 
behavior: 

DELSIG  -  CE  x  DELEPS 

(c)  Calculate  TAU: 

TAU  «  SIG  *  DELSIG 

(d)  With  TAU  as  the  state  of  stress,  determine  the  value  of  the  yield 
function  F. 

(e)  If  F{TAU)  <  0,  elastic  behavior  assumption  correct  (loading  elastically, 
neutral  loldlng,  or  unloading).  Hence,  TAU  Is  the  correct  stress  value 
and  the  calculations  for  this  algorithm  are  complete  (RETURN).  If 
F(TAU)  >  0,  steps  (f),  (g),  and  (h)  must  be  executed. 

(f)  If  the  previous  state  of  stress  was  plastic  (as  indicated  by  a  flag),  set 
RATIO  «  0  and  go  to  step  (g).  Otherwise,  there  Is  a  transition  from 
elastic  to  plastic  and  RATIO,  the  portion  of  DELSIG  which  Is  elastic, 
must  be  determined.  The  variable  RATIO  Is  determined  from  the  equation 

F  {SIG  ♦  (RATIO  x  DELSIG))  -  0 

since  at  the  stress  SIG  ♦  (RATIO  x  DELSIG)  the  yield  function  F  -  0  and 

yielding  is  initiated. 

(g)  Redefine  TAU  as  the  stress  at  start  of  yield 

TAU  -  SIG  ♦  (RATIO  x  DELSIG) 
and  calculate  the  elastic-plastic  strain  increment 
OEPS  -  (1  -  RATIO)  x  DELEPS 

(h)  To  obtain  the  final  stresses,  which  include  the  effect  of  the  complete 
strain  increment  DELEPS,  the  stresses  corresponding  to  the  elastic- 
plastic  strain  increment  DEPS  must  be  added  to  TAU.  Hence,  TAU  Is  the 

current  state  of  stress  {c)n+1 

TAU  -  TAU  ♦  (CEP  x  OEPS) 

and  t*e  algorithm  calculations  are  complete  (RETURN). 


t 

I 


were  determined  in  STRES3  based  on  the 


Incremental  stresses,  (Ac), 
previously  calculated  strain  increments  (step  b  in  Table  4.6).  These 
stresses  were  computed  at  elemental  quadrature  points  as  was  done  for 
the  velocity  strains.  These  incremental  stress  calculations  were 
performed  under  the  assumption  that  the  strain  increments  were  entirely 
elastic.  For  this  assumption,  the  incremental  stresses  were  computed 
using  elastic  stress-strain  relations.  Equations  4.19  and  4.20  show  the 
elastic  stress-strain  matrix  used  in  the  determination  of  the 
incremental  stresses  and  the  corresponding  multiplication  process  which 
was  executed.  The  incremental  strain  in  the  z -direction,  c2,  was 
computed  according  to  Equation  4.21. 


(1  -  v  ) 


(4.19) 


{Aa}  =  C  X  {Ac} 


(4.20) 


where. 


(Aa}T 


{Ac}1 


[a<jx  Acy  at xy],  and 

AAA 

[Acx  Aey  Ayxy:. 


Acz  «  (Acx  +  Acy) 


(4.211 


The  total  stresses  for  the  current  time  step,  {o)n+1,  wen 
calculated  by  updating  the  total  stresses  from  the  previous  time  step, 
f 3 } n ,  as  shown  in  Equation  4.22. 
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{o}n+*  =  {a}n  +  (Ao)  (4.22) 

These  updated  total  stresses,  {a}n+*,  were  then  implemented  into  the 
corresponding  yield  function  (F)  in  order  to  determine  if 
elastic-plastic  behavior  had  occurred  during  the  time  step.  The  yield 
function  used  in  STRES3  for  the  plane  stress  material  law  is  shown  in 
Equation  4.23. 

F  =  (°x  +  °y  "  ax  ay  +  3Txy}i  (4‘23) 

This  function  was  based  on  the  von  Mises  yield  criterion  (11,  12,  15) 
and  was  evaluated  in  subroutine  STREQU  using  values  obtained  from 
STRES3.  The  value  for  F  determined  in  Equation  4.23  was  compared  to  the 
yield  stress  value  (a^ld)  obtained  from  the  previous  time  step.  This 
yield  stress  value  was  initially  established  such  that  where 
a1  =  the  input  stress  value  for  the  first  point  (ej.o^)  on  the 
monotonically  increasing  piecewise  linear  stress-strain  curve.  This 
value  for  a^d  remained  unchanged  for  a  particular  elemental  quadrature 
point  until  yielding  was  initiated.  The  yield  stress  value  was  then 
redetermined  such  that  a^d  *  F(an),  except  for  the  case  of 
elastic-perfectly  plastic  material  behavior,  where  remains  equal  to 
<jy  The  yield  value,  o^ld,  was  set  equal  to  F(on)  in  order  to  record 
the  previous  state  of  stress  or  stress  history  for  use  in  comparison 
with  the  current  state  of  stress  as  determined  by  Equation  4.23.  In 
comparing  F  to  <7^d»  if  F  <  c”^d,  then  the  assumption  of  elastic 
behavior  (loading  elastically,  neutral  loading,  or  unloading)  was 
correct.  Therefore,  the  stresses  computed  in  Equation  4.22  were  the 


1 
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correct  total  stresses  for  the  current  time  step  (step  e  in  Table  4.6). 
However,  if  F  >  then  elastic-plastic  behavior  has  occurred  and  the 
stresses  obtained  from  Equation  4.22  were  incorrect  and  were 
recalculated  according  to  the  elastic-plastic  stress-strain  relations  in 
the  following  paragraphs. 

Two  additional  calculations  were  necessary  for  elastic-plastic 
material  behavior.  The  first  calculation  was  the  determination  of  the 
plastic  modulus  (modulus  of  the  next  segment  on  the  piecewise  linear 


stress-strain  curve)  and  the  next  yield  point.  The  plastic  modulus  (Ey) 
was  determined  in  subroutine  STRMOD  using  Equation  4.24  and  the  yield 


stress  yalue  was  set  equal  to  a j  such  that  >  0^^. 


c  _  qI  '  °I-1 
T  £I  "  eI-l 


(4.24) 


where. 


aI,cI 


'i-r€i-i  “ 


plastic  modulus, 

stress  and  strain  values  for  the  next  yield  point  on 
the  linear  stress-strain  curve,  and 

stress  and  strain  values  for  the  current  yield 
point. 


However,  if  the  stress-strain  point  Uj.cjj)  did  not  exist,  perfect 
plasticity  was  assumed  (Ey  3  0)  and  the  yield  stress  value  remained 
equal  to  Oj  j.  As  a  note,  the  initial  von  Mises  yield  surface  that  was 
established  for  the  initial  yield  stress  was  simply  expanded  for  the 
following  yield  stresses  using  an  isotropic  hardening  rule. 


I 


i 


The  second  calculation  performed  was  the  evaluation  of 
Equation  4.23  using  the  stresses  from  the  previous  time  step.  This 
value  F(an)  was  then  compared  with  o^d  using  Equation  4.25. 


F(on)  -  o" 


SHIST  = 


(4.25) 


Equation  4.25  was  evaluated  in  order  to  determine  if  the  previous  state 
of  stress  was  elastic  or  plastic.  If  SHIST  <  0.001,  then  the  previous 
state  of  stress  was  plastic  and  the  incremental  strains  calculated  in 
Section  4.2.3  were  entirely  plastic.  Therefore,  the  incremental 
stresses  calculated  using  Equation  4.20  were  Invalid.  These  incremental 
stresses  were  then  recalculated  using  an  elastic-plastic  stress-strain 

CD 

matrix  C  .  However,  if  SHIST  >  0.001,  then  the  previous  stress  state 

A 

was  elastic  and  the  portion  (RATIO)  of  the  Incremental  stresses,  (Ao), 
which  were  elastic  was  determined  us'  g  Equation  4.26.  This  equation 
was  obtained  based  on  the  expression  shown  In  Equation  4.27. 


RATIO 


where , 


t  *  A1  *  4fF(a»)a  «  (<>;,„»  -  r(»")»)j 

2F(Ac)2 


-2Ao  on  -  2Ao  cn  +  Ao  cn  ♦  Ac  on  -  6At  rn 
x  x  y  y  yx  xy  xyxy 


(4.26) 


F[on  +  (RATIO  x  Ac)]  -  cj]d  •  0 


(4.27) 


Equation  4.26  was  verified  by  comparison  with  the  solution  to  4.27.  The 
subroutine  STRYIE  was  used  to  calculate  RATIO.  The  elastic  portion  of 


the  incremental  stresses  was  calculated  using  RATIO  multiplied  by  the 
previously  calculated  incremental  stresses.  These  elastic  incremental 
stresses  were  then  added  to  the  stresses  from  the  previous  time  step  as 
shown  in  Equation  4.28. 


(c}n  =  {a}n  +  RATIO  x  (Ac) 


(4.28) 


Equation  4.28  was  executed  in  order  to  update  the  old  stresses  to 

include  the  elastic  portion  of  the  incremental  stresses.  These  updated 

stresses  were  then  implemented  into  the  yield  function  in  order  to 

update  the  stress  state  to  include  the  elastic  stress  Increments, 

°yld  =  F(°n)-  This  yield  function  value  was  used  later  in  the 

rp 

determination  of  the  elastic-plastic  stress-strain  matrix  C  .  The 
plastic  strain  increments  (Aexp,  4eyp,  and  Arxyp)  and  the  elastic  strain 

A 

increment  in  the  z-direction  (Aez)  were  calculated  using  RATIO  as  shown 
in  Equations  4.29  and  4.30. 


<v 


=  (1  -  RATIO)  ( Ac) 


where. 


{Aep} 


plastic  strain  increments,  and 
incremental  strains  from  Section  4.2.3. 


(4.29) 


Aez  *  RATIO  x  Ae2 


(4.30) 


The  plastic  stress  increments  were  calculated  using  elastic-plastic 

stress-strain  relations  developed  in  accordance  with  the  yield  function 

(F),  isotropic  hardening  and  an  associated  flow  rule.  These  elastic- 

EP 

plastic  stress-strain  relations  are  incorporated  into  the  C  matrix 
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shown  in  Table  4.7.  This  matrix  was  developed  for  incremental 

plasticity  analyses  using  the  yield  function  from  Equation  4.23  and  an 

EP 

associated  flow  rule.  The  derivation  of  C  for  the  plane  stress 
material  law  was  similar  to  that  shown  in  Appendix  C  for  the  biaxial 
elastic-plastic  plane  strain  or  axisymmetric  material  law.  The 
incremental  plastic  stresses  were  determined  using  Equation  4.31.  The 
incremental  strain,  both  the  elastic  and  plastic  portions,  in  the 
z-direction,  Aez,  was  calculated  using  Equation  4.32.  However, 
Equation  4.32  is  not  correct.  The  correct  formula  for  determining 
is  shown  in  Equation  4.33. 


Uop} 


CEP  x  Uep} 


(4.31) 


where. 


<v 


{Aep> 


incremental  plastic  stresses,  (Aoxp  &o^p  ATxyp), 


elastic-plastic  stress-strain  matrix  shown  in 
Table  4.7,  and 

AAA 

plastic  strain  increments,  (Aexp  Aeyp  Ayxyp)  which 

were  either  the  strains  calculated  in  Section  4.2.3 
(previous  state  of  stress  plastic)  or  the  strains 
calculated  using  Equation  4.29. 
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— ~6  Sl2Sf  4yxv  )  (1  -  RATIO) 
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where, 


Ae  ,  Ae  .  Ay  =  incremental  strains  calculated  in 
y  ^  Section  4.2.3, 


3  1/  1 

ipr-j  < — n 

yld  1  +  yn 


-)  ,  and 


=  value  from  Equation  4.21  if  RATIO  *  0  or 

value  from  Equation  4.30  (Ae^)  if  RATIO  f  0. 

.  ii  .  {  ^  '  8  S11S33>  s-  +  (T^  ~  S  S22S33>  ^  + 

Z  ,  1-v  .  B  s?  ■  x  ~  y 


/i-v  .  B  sf-, 
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•6  S1?S33 

- -  }  (1  -  RATIO) 

fl-v  -  B  S*  ** 

(ir^j  33) 


(4.33) 


The  only  differences  between  the  two  equations  exist  in  the  expressions 
containing  Poisson's  ratio.  Equation  4.32  contains  the  terms  -p^  and 
-pj  while  the  corresponding  terms  in  Equation  4.33  are  pp  and  pp. 


The  terms  inside  the  brackets  on  the  righthand  side  of  Equation  4.33 

A 

were  obtained  by  substituting  Ac^p  *  0  in  the  elastic-plastic 
stress-strain  relations  for  the  plane  strain  or  axisymnetric  material 
law  and  solving  for  Aezp. 

The  total  stresses  for  the  current  time  step  (a}n+*  were  determined 
using  the  plastic  increments  of  stress  obtained  from  the  evaluation  of 
Equation  4.31.  The  total  stresses  were  calculated  using  Equation  4.34. 


(a}n+1  =  {o}n  +  {Ao  } 

P 


(4.34) 


where, 


V 

s 

5 


-  ‘ -.•'o' o' oV ' 


stresses  for  the  current  time  step. 


either  the  stresses  from  the  previous  time  step  if 
RATIO  3  0  or  the  stresses  calculated  in 
Equation  4.28  if  RATIO  f  0,  and 

plastic  increments  of  stress  determined  using 
Equation  4.31. 

A  value  for  the  yield  function  (F)  was  then  calculated  with  these  new 
stresses,  F(on+1).  This  value  F(an+1)  was  then  compared  with  the  value 
for  the  next  yield  point,  Oj,  obtained  from  subroutine  STRMOO.  This  was 
executed  in  order  to  determine  If  the  current  state  of  stress  exceeded 
the  failure  criterion  for  the  next  yield  point  Oj.  If  F(an'*’*)  <  Oj, 
then  the  stresses  calculated  using  Equation  4.34  were  correct  and 
ayld  *  F(°n+1)-  However,  if  F(an+1)  >  Oj,  then  F(on)  »  and 

°yld  *  °i»  and  th€  calculation  process  was  reinitiated  beginning  with 
the  determination  of  the  next  value  for  the  plastic  modulus 
(Equation  4.24)  and  the  next  yield  point  value.  One  other  possibility 
existed  for  the  case  of  perfect  plasticity.  For  this  case,  a  reduction 
factor  RED  was  determined  using  Equation  4.35  which  was  then  multiplied 
by  the  stresses  obtained  from  Equation  4.34  as  shown  in  Equation  4.36. 
After  these  two  equations  were  executed.  Equation  4.37  was  evaluated. 


(Ac  _} 


RED 


Jyld 


F(on+1) 


(4.35) 


where, 


RED 


factor  used  to  reduce  the  stresses  calculated  in 
equation  4.34  so  that  the  failure  criterion, 

F(on+1)  <  oJld,  was  met. 


A 


{a}n+1  =  RED  x  {a}n+1 


SHIST2  = 


F(an+1)  -  a" 


(4.36) 


(4.37) 


where. 


(a}n+1 


stresses  obtained  from  Equation  4.36. 


If  the  value  for  SHIST2  was  less  than  or  equal  to  0.001  the  stresses 
calculated  using  Equation  4.36  were  the  correct  total  stresses  for  the 
current  time  step.  But,  if  SHIST2  was  greater  than  0.001, 
Equation  4.35,  4.36,  and  4.37  were  reexecuted  until  SHIST2  was  found  to 
be  less  than  or  equal  to  0.001. 

Upon  the  determination  of  the  total  stresses  for  the  current  time 
step,  the  total  strains  were  updated  to  include  the  incremental  strains 
as  shown  in  Equation  4.38.  These  new  strains  and  the  new  stresses  were 
then  stored  in  the  appropriate  location  of  the  strs  array.  The  value 


for  c"]d  was  also  stored  in  the  strs  array. 


(c}n+1  «  (e}n  +  Ue) 


(4.38) 


where , 


(4e)  »  strains  obtained  from  Section  4.2.3  (4c  4e  and  4y  ) 
and  the  strain  (4ez)  computed  in  STRES3.  y  y 

The  entire  formulation  contained  in  STRES3  for  the  biaxial 
elastic-plastic  plane  stress  material  law  was  found  to  be  correct, 
except  for  Equation  4.32  which  was  used  to  determine  the  total  strain 
increment  in  the  z-directlon  (€z).  However,  it  should  be  mentioned  that 
there  should  be  more  Information  in  the  SAMS0N2  Users  Manual  (3)  with 


oVv-AV-j.'- 


WPS 


{ 
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regard  to  the  types  of  materials  that  can  be  correctly  modeled  using  the 
STRES3  subroutine.  STRES3  should  only  be  used  for  problems  that  involve 
materials  which  exhibit  ductile  or  semiductile  stress-strain  material 
behavior.  In  fact,  the  von  Mises  failure  criterion  is  only  good  for  a 
yielding  mode  of  failure  which  is  exhibited  in  ductile  or  semiductile 
metals. 

4. 2.4.3  Stress  Evaluations  Using  the  Biaxial  Elastic-Plastic 
Plane  Strain  or  Axisymmetric  Material  Law  (STRES3) 

The  formulation  for  the  bia/ial  elastic-plastic  plane  strain  or 
axisymmetric  material  law  was  also  contained  in  subroutine  STRES3. 
Therefore,  the  identical  solution  scheme  (Table  4.6)  was  used  for  the 
plane  strain  or  axisymmetric  material  law  as  was  used  for  the  plane 
strain  material  law  with  a  few  changes  in  the  equations  used.  In  this 
section,  due  to  this  similarity  in  the  solution  scheme,  only  the 
equations  which  were  different  are  displayed.  The  discussion  of  the 
solution  algorithm  will  not  be  presented  however. 

The  first  difference  between  the  two  material  laws  was  in  the 
elastic  stress-strain  matrix,  C^.  Equations  4.39  and  4.40  show  the 
matrix  and  the  corresponding  multiplication  process  associated  with  C^. 

,E  _  E 

c  -  TTnTTT^T 


1 

I 


1-v 


V 

0  0 


1“V  0  V 

l-2v 


0 
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(4.39) 


{Aa}  =  C  x  {Ac} 


(4.40) 


where, 

{Ao)T  =  [Aa  Aa  At  Ac  ]  ,  and 

*  y  z 

{Ae}^"  *  [Ac  Ac  Ay  Ac  ]. 

L  x  y  'xy  zJ* 

As  illustrated  in  these  two  equations,  one  additional  stress  term  (Ac  ) 
and  one  additional  strain  term  (A c^)  were  present.  However,  the  Ac^ 
term  equalled  zero  for  plane  strain  analyses.  It  should  be  noted  that 
the  Equations  4.21,  4.30,  4.32,  and  4.33  that  were  used  to  determine  the 
strain  increment  in  the  z-direction  (a t^)  for  the  plane  stress  material 

law  do  not  apply  to  the  plane  strain  and  axisymmetric  material  law.  The 

A ez  term  for  axisymmetric  analyses  was  calculated  in  subroutine  V8FRCN. 

The  second  difference  between  the  two  material  laws  was  the  yield 
function  that  was  used.  The  yield  function  for  the  plane  strain  or 
axisymmetric  material  law  is  shown  in  Equation  4.41.  Three  additional 
terms  are  present  in  Equation  4.41  compared  to  Equation  4.23.  These 

A 

terms  are  all  a  function  of  cz  which  equalled  zero  for  plane  stress. 

F  =  (a2  +  a2  ♦  a2  -  a  a  -  a  a  -  a  a,  +  3t2  (4.41) 

x  y  z  xy  xz  yz  xy' 

The  third  difference  that  existed  between  the  two  material  laws 
involved  the  equation  for  determining  RATIO.  The  value  for  t  in 
Equation  4.26  needed  to  be  updated  in  order  to  include  the  o2  terms. 
Equation  4.42  shows  the  formula  that  was  used  to  determine  t  for  the 
plane  strain  or  axisymmetric  material  law. 
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(4.42) 


The  value  for  t  determined  in  Equation  4.42  was  substituted  into 
Equation  4.26  in  order  to  determine  which  portions  (RATIO)  of  the  stress 
and  strain  increments  were  elastic. 


The  final  difference  between  the  two  material  laws  existed  in  the 

.EP 


elastic-plastic  stress-strain  matrix  C 


The  derivation  involved  in 


EP  EP 

the  determination  of  C  and  the  final  C  matrix  for  the  plane  strain 


or  axisymmetric  material  law  is  shown  in  Appendix  C.  The  incremental 


stresses  and  strains  corresponding  to  Equation  4.31  are  {Acp}  =  ^&cXp 

AM  A  A  T  AAA  A 

Ac  At  Ac  1  and  (Ac  }  a  \ Ae  Ae  Ay  Ae  1. 

yP  xyp  °°zpJ  p  -  xp  0  yp  Txyp  otzpJ 


The  remainder  of  the  calculations  performed  in  STRES3  were 
identical  for  the  two  material  laws.  But,  there  were  four  stresses  and 
four  strains  involved  in  the  calculations  for  the  plane  strain  or 
axisymmetric  material  law  and  only  three  stresses  and  three  strains  for 
the  plane  stress  material  law. 

The  formulation  for  the  plane  strain  or  axisymmetric  material  law 
was  found  to  be  correct  through  comparisons  with  References  12,  14,  and 
15.  No  errors  were  found  in  the  entire  plane  strain  or  axisymmetric 
formulation. 


4. 2. 4. 4  Stress  Evaluations  Using  the  Plane  Strain  or  Axisymmetric 
AFWL  "Engineet  ,nq"  Material  Law  (STRES9) 

This  section  contains  a  brief  discussion  on  the  AFWL  "engineering" 
model.  In  addition,  some  of  the  principal  equations  that  are  used  in 
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this  model  are  displayed.  The  AFWL  "engineering"  material  law  was 

developed  primarily  for  modeling  soils.  It  is  defined  by  a  piecewise 

linear  hydrostat  (plot  of  hydrostatic  pressure  ([a  +  a  +  a  ]/3)  versus 

x  y  2 

volumetric  strain)  with  corresponding  bulk  moduli  for  both  loading  and 
unloading/reloading  («L  and  Ky)  (see  Figure  4.3). 


Lockup 


Unloading  /  Reloading 


Load i no 


Tension  Cutoff 


Volumetric  Strain 


Figure  4.3  An  Example  of  a  Piecewise  Linear  Hydrostat  Curve  Which  is 
Used  to  Define  the  AFWL  "Engineering"  Material  Law. 

It  is  also  defined  by  a  yield  (failure)  surface  which  plots  hydrostatic 
pressure  versus  ^  (Drucker-Prager  yield  criterion)  where  ^  equals 
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the  square  root  of  the  second  invariant  of  the  deviatoric  stress  tensor. 
A  nonassociated  flow  rule  is  employed  in  the  AFWL  "engineering"  model. 
It  should  be  mentioned  that  the  AFWL  "engineering"  model  can  only  be 
used  to  model  strain  hardening  materials.  Also,  the  AFWL  model  was 
formulated  only  for  plane  strain  or  axisymmetric  analyses.  Hence,  the 
stresses,  (a),  and  strains,  {cl,  that  were  used  in  the  calculations  are 


as  follows: 

A  V  AAA  A  A  y  AAA  A 

{ o }  *  To  o  t  o  1  and  (e)  =  [o  a  x  o  ]. 

-  x  y  xy  zJ  L  x  y  xy  zJ 

Some  preliminary  operations  were  performed  prior  to  the 
determination  of  the  stress  values  at  the  elemental  integration  points. 
The  first  operation  was  the  determination  of  the  total  strains  for  the 

A  _  ,  1 

current  time  step,  { e >  1  (see  Equation  4.38).  The  next  operation  was 

the  determination  of  the  four  history  parameters  which  are  described 
below: 


pvolp 

pevol 

d4 

d6 


hydrostatic  pressure  for  the  previous  time  step, 

Pn, 

volumetric  strain  for  the  previous  time  step,  e"Ql, 

minimum  volumetric  strain  which  has  been  computed 
throughout  the  solution,  and 

volumetric  strain  at  which  the  tension  cut-off  was 
entered. 


.*> 


» 


The  final  operation  was  the  calculation  of  the  volumetric  strain  for  the 

*  _  i  1 

current  time  step  (cy0{)  as  shown  in  Equation  4.43. 


(4.43) 
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The  value  which  was  obtained  from  Equation  4.43  was  compared  with  the 
history  parameter  d4  in  order  to  determine  which  solution  scheme  was  to 
be  executed.  There  are  two  possible  solution  schemes  that  are  used  in 
STRES9.  The  first  is  for  a  loading  case  where  <_  d4  and  the  second 
scheme  is  for  an  unloading/reloading  case  >  d4.  These  tro  schemes 
are  discussed  separately  in  the  following  paragraphs. 

If  the  current  value  for  the  volumetric  strain  was  less  than  or 
equal  to  the  value  for  the  history  parameter  d4,  then  the  solution  was 
determined  to  be  in  a  state  of  loading.  Upon  this  determination,  a  bulk 
modulus  (K  *  K^)  value  was  determined  from  the  hydrostat  curve  based  on 

A  _  i  1  A  _  ,  1 

the  value  for  e^.  Also,  the  value  for  d4  was  set  equal  to  e"  j  and  a 

value  for  v  (Poisson's  ratio)  was  chosen  from  the  input  data  based  on 

evol‘  The  normal  and  shear  stresses  for  the  current  time  step  (o"+1  and 

t^1)  were  then  calculated  for  particular  integration  points  according 
xy 

to  Equations  4.44,  4.45  and  4.46. 

“n+i  _  *n 

+  26  [«, .  ( xaU]  +  4P  («.««) 

where , 

2G  .  3K. 

evol  "  pevo1*  and 

aP  *  incremental  change  in  the  hydrostatic  pressure  between 

the  current  and  previous  time  steps  as  was  calculated 

using  Equation  4.46. 
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These  three  equations  were  found  to  be  correct  when  compared  to 


information  that  was  obtained  from  AFWL  personnel.  However,  a  possible 


error  was  discovered.  The  stress  in  the  z-direction  (c  )  was  set  equal 


to  0.0  for  elements  having  a  thickness  not  equal  to  1.0.  This  appears 


to  be  a  plane  stress  correction.  The  current  value  for  the  hydrostatic 


pressure  was  calculated  using  Equation  4.47. 


>n+l  =  p"  + 


(4.47) 


where. 


time  history  parameter  pvolp. 


This  value  for  the  current  hydrostatic  pressure  was  then  compared  to  the 


limiting  mean  normal  stress  value  P„..M  (the  tension  cut-off  value).  Two 

max 


possible  cases  were  obtained  based  on  this  comparison.  The  first  case 


was  for  Pn+1  -  <0.  This  case  showed  that  the  loading  had  occurred 

max 


outside  of  the  tension  cut-off  region  of  the  hydrostat  curve.  For  this 

case  (Pn+1  -  P„v  <  0),  a  value  for  the  failure  function  (XJ)  was 
ma  x 

computed  by  substituting  the  stress  values  that  were  calculated  using 


Equations  4.44  and  4.45  into  Equation  4.48. 


z  '  '  y 


on+V  ♦  6(Tn"V] 
z  '  '  xy  '  J 


(4.48) 


The  value  XJ  was  then  compared  with  a  value  YJ  which  was  obtained  from 
the  Drucker-Prager  (17)  yield  (failure)  surface  (plot  of  ^  versus 
hydrostatic  pressure).  The  value  for  YJ  was  s^t  equal  to  the  value  for 


which  was  obtained  from  the  yield  surface  corresponding  to  the  value 


for  Pn+1.  The  comparison  between  YJ  and  XJ  was  performed  in  order  to 
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determine  if  yielding  had  occurred.  If  XJ  <  YJ,  then  the  stresses  which 
were  previously  computed  in  Equations  4.44  and  4.45  were  correct. 
However,  if  XJ  >  YJ,  then  the  nonassociated  flow  rule  that  was 
formulated  in  STRES9  was  applied  and  the  stresses  were  recalculated  by 
substituting  the  previous  stress  values  from  Equations  4.44  and  4.45 
into  Equation  4.49. 


.  pn+1  *  com  <;j+1  -  pntli 

;n+1  •  com  ;n+1 
xy  xy 


(4.49) 


where. 


i  YJ 
coni  =  y-j. 

The  stress  values  that  were  obtained  from  Equation  4.49  were  the  correct 

values  for  the  case  where  yielding  had  occurred.  The  second  case  which 

was  obtained  from  the  comparison  between  Pn+1  and  Pmax  was  for 

Pn+1  -  P  .  >0.  This  comparison  shows  that  the  current  state  of  stress 

occurs  beyond  the  tension  cut-off  region.  Therefore,  the  values  for  the 

stresses  that  were  computed  in  Equations  4.44  and  4.45  exceed  the 

tension  cut-off  value  and  must  be  corrected.  These  values  were  changed 

such  that  all  of  the  normal  stress  values  (ax+1,  °y+*»  and  °2+1)  and  the 

value  for  Pn+*  were  set  equal  to  the  limiting  mean  normal  stress  value 

P  .  The  shear  stress  value  was  set  equal  to  0.0.  The  time  history 
max 

parameter  d6  was  also  modified  using  Equation  4.50. 


*  el,  ♦  (P, 


vol  '  max 


Pn)  /  K 


(4.50) 


where. 


d6  *  value  of  the  volumetric  strain  when  the  solution  entered 
the  tension  cut-off  region. 
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If  the  current  value  for  the  volumetric  strain  ^  was  found  to  be 
greater  than  the  value  for  d4,  then  the  solution  was  determined  to  be  in 
an  unloading/reloading  state.  The  value  for  the  bulk  modulus  (X  =  K^) 
and  the  value  for  Poisson's  ratio  (v)  were  determined,  as  in  the  loading 
case,  based  on  the  value  for  e"0|  as  long  as  pvolp  <  pmax.  These  values 
for  the  bulk  modulus  and  Poisson's  ratio  were  then  substituted  into 
Equations  4.44,  4.45  and  4.46,  as  was  done  for  the  loading  case,  and  the 
remainder  of  the  solution  scheme  was  identical  to  that  for  the  loading 
case.  However,  if  pvolp  >  P  ,  then  another  comparison  was  required 

**■  TTVo  X 

between  and  d6.  If  was  greater  than  or  equal  to  d6,  then  the 
normal  stress  values  remained  equal  to  P_„„  and  the  shear  stress  values 
equal  0.0.  For  the  other  case  where  <  d6,  the  stresses  were 
recalculated  following  the  same  procedure  mentioned  earlier  in  this 
paragraph  for  the  loading/reloading  case  where  pvolp  <  P  .  However, 
the  formulation  in  subroutine  STRES9  was  found  to  be  incorrect  for  this 
case  where  pvolp  >  Pmax  and  <  d6  when  compared  to  information  which 
was  obtained  from  personnel  at  the  AFWL.  The  total  increment  of 
volumetric  strain  cannot  be  used  in  order  to  determine  the  stresses  for 
a  case  where  reloading  has  occurred  beyond  the  tension  cut-off.  The 
proper  formulation  requires  the  determination  of  a  factor  Fa  which  is 
equal  to  the  percent  of  the  total  strain  Increment  which  is  actually 
reloading  out  of  the  tension  cut-off.  This  factor  F#  is  computed  using 
Equation  4.51  and  needs  to  be  implemented  Into  Equations  4.44,  4.45,  and 
4.46  in  order  to  calculate  the  correct  values  for  the  stresses. 


.'•W ,vv.\*v 


112 


F 

a 


n+1 

£vol 
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Evol 


(4.51) 


Two  additional  operations  were  performed  after  the  correct  stress 
values  were  obtained  for  either  the  loading  case  or  the  unloading/ 
reloading  case.  The  first  operation  which  was  executed  was  the  storage 
of  the  current  stress  and  strain  values,  which  were  calculated  at 
particular  integration  points,  in  the  strs  array.  The  second  operation 
was  the  storage  of  the  values  for  the  four  history  parameters  in  the 
strs  array. 

The  formulation  for  the  AFWL  material  model  (STRES9)  was  verified 
through  comparisons  with  information  that  was  obtained  from  AFWL  .  _ 
personnel  and  two  errors  were  found.  Also,  it  was  found  that  a  segment 
within  the  STRES9  subroutine  was  never  used  in  any  analysis.  This 
segment  was  associated  with  a  plane  stress  material  law  and  appears  to 
have  been  eliminated  from  use  in  the  AFWL  model. 


4. 2. 4. 5  Evaluation  of  the  Internal  Hodal  Forces 

This  section  contains  the  development  of  the  equations  used  in 

SAMS0N2  for  the  determination  of  the  internal  nodal  forces  'W- 

These  forces  are  produced  as  a  result  of  the  stresses  and  strains  which 

develop  within  an  element  (see  Table  4.4  on  page  83).  Therefore,  the 

internal  nodal  forces  (F.  .}  are  a  function  of  the  amount  of  stress  and 

int 

strain  exhibited  within  an  element. 

Equation  4.52  is  the  basic  equation  used  in  the  determination  of 
the  internal  nodal  forces.  However,  no  stiffness  matrix  is  ever 


determined  in  the  execution  of  SAMS0N2  as  was  mentioned  in  Section  4.1. 
Therefore,  the  terms  on  the  righthand  side  of  Equation  4.52  need  to  be 
changed  into  terms  that  are  consistent  with  the  SAMS0N2  formulation. 

<Fint>  =  [K]{U)  (4.52) 
where. 


(F.  J 
int 

- 

internal  nodal  force 

vector , 

[K] 

3 

structural  stiffness 

matrix,  and 

(U) 

- 

vector  containing  the  nodal  displacements. 

The  stiffness  matrix,  [K],  can  be  expressed  in  terms  of  the  strain- 
displacement  and  stress-strain  matrices  ([B]  and  [C])  that  were 
calculated  in  the  previous  sections  as  given  in  Equation  4.53. 

[K]  *  /  [B]T  [C]  m  dV  (4.53) 

v 

where , 

[K]  *  elemental  stiffness  matrix  calculated  at  a  particular 
integration  point, 

[B]  *  [B]/IJl  ,  where  [B]  is  the  strain-displacement  matrix 

shown  in  Table  4.5  of  Section  4.2,3,  and 

[C]  *  stress-strain  matrix  for  a  particular  material  model. 

The  next  step  is  to  substitute  this  expression  for  [K]  into 
Equation  4.52  which  results  in  Equation  4.54. 

(F.  .}  *  /  [B]T  [C]  CB1  dV  (U>  (4.54) 

I  nt  y 

But,  this  equation  for  {F}  can  be  simplified  using  the  expressions  in 
the  following  equation: 
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U)  =  [B]  {U> 

r  ,  «.55) 

{0}  *  [C]  {c} 

The  equation  resulting  from  the  substitution  of  these  two  expressions 
into  Equation  4.54  is  shown  in  Equation  4.56.  However,  Equation  4.56 
still  needs  to  be  altered  in  order  to  incorporate  the  Gauss-Legendre 
numerical  integration  scheme. 

{F}  =  /  [B]T  (a)  dV  (4.56) 

v 


Equations  4.57  and  4.58  are  the  final  expressions  that  were  formulated 
into  the  SAMS0N2  code  in  subroutine  V8FRCN  for  the  evaluation  of  the 
internal  nodal  forces  for  plane  and  axisymmetric  analyses,  respectively. 

n  n 


n+1 


t  i  z  h.  h.  [BU-.n,)]  (oU-.nJ} 
i*l  j*l  1  J  1  J  1  J 


n+1 


(4.57) 


where. 


<Fint> 


n+1 


t 


Vhj 

CB({1.nj)]T 


«  internal  forces  applied  to  the  nodes  of  an  8NQ 
element  for  a  plane  analysis, 

“  Cpxl  Fyl  Fx2  Fy2  *•*  Fx7  Fy7  Fx8  Fy8^’  where 
Fxl*  Fyl  *  lnternal  nodal  forces  in  the  x  and  y 

direction  for  node  I, 

*  thickness  of  the  element, 

■  order  of  Integration  selected  for  the  analysis, 

*  Gauss-Legendre  weight  coefficients, 

»  transpose  of  the  stress-strain  matrix  shown  in 
Table  4.5  which  was  eval"ated  at  a  particular 
e'emental  quadrature  point  with  the  coordinates 
U^.hj)  being  assigned  values  corresponding  to  the 

Gauss-Legendre  abscissae  coefficients,  and 
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{a(£. ,r  ■) }n+*  =  total  stresses  for  the  current  time  step,  including 
J  the  viscous  stresses  due  to  damping,  which  were 

evaluated  at  elemental  quadrature  points  using  an 
appropriate  material  law  and  were  then  transformed 
into  the  global  (x,y)  coordinate  system. 


{Fint} 


n+1 


^  ^  h.hj  r(Ci,nj)[B(Ci,nj)]T{a(Ci,nj)}n+1 


(4.58) 


where, 

(F.  }n+1  *  internal  forces  applied  to  the  nodes  of  an  8NQ 

int  element  for  an  axisymmetric  analysis,  and 


r(£.,n.)  =  distance  from  the  axis  of  symmetry  (y-axis)  to  a 

1  J  particular  quadrature  point  (5,,n.)  obtained  by  using 

'  J 

Equation  4.1  in  Section  4.2.1  and  solving  only  for  x. 


Equations  4.57  and  4.58  were  developed  by  replacing  the  integral  in 
Equation  4.56  by  the  appropriate  terms  used  in  the  Gauss-legendre 


integration  scheme.  In  addition,  the  expressions  in  Equations  4.59  and 


4.60  were  used  to  simplify  Equation  4.56. 


fm 


3 


dV  *  r°  | J |  d£ 

dn 

(4.59) 

where. 

d  V 

differential  volume  used  in  the  integration. 

r* 

s 

t  for  plane  analyses  (a  *  0),  or 

r*{ ci  ,n . )  for  axisymmetric  analysis  (a  « 

1),  and 

d£  dn  « 

differential  area  used  in  the  numerical 
scheme. 

integration 

[8]T  >  Ul  [8]' 

T 

(4.60) 

It  should  be  noted  that  the  total  internal  nodal  forces  for  a  particular 
node  included  the  contributions  from  all  of  the  elements  to  which  it  was 
associated. 
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The  stresses  that  were  used  in  Equations  4.57  and  4.58  were  not 
equal  to  the  stresses  that  were  calculated  in  the  appropriate  material 
law  subroutines.  The  stresses  which  were  calculated  using  the  material 
laws  were  first  adjusted  to  include  the  viscous  stresses  that  were 
calculated  in  order  to  simulate  stiffness  proportional  damping.  These 
adjusted  stresses  were  then  rotated  from  the  corotational  (x,y) 
coordinate  system  to  the  global  (x,y)  coordinate  system  by  using  a 
suitable  transformation  matrix.  Both  of  these  adjustments  to  the  stress 
values  that  were  obtained  from  the  material  law  subroutines  are 
discussed  in  the  following  paragraph.  One  other  observation  that  was 
made  regarding  only  Equation  4.58  was  that  the  factor  2*  was  not 
included.  However,  this  factor  is  canceled  later  in  obtaining  the  final 
solution  to  the  equation  of  motion. 

The  stresses  (a}n+1  that  were  calculated  previously  in  the  material 
law  subroutines  were  modified  twice  prior  to  their  use  in  Equations  4.57 
and  4.58.  The  first  modification  was  the  addition  of  the  viscous 
stresses  {a}vis  which  were  calculated  in  SAMS0N2  In  order  to  simulate 
stiffness  proportional  damping  through  the  use  of  linear  artificial 
viscosity.  The  reason  for  the  use  of  linear  artificial  viscosity  is 
because  a  stiffness  matrix  is  never  created  in  the  SAMS0N2  explicit 
code.  Linear  artificial  viscosity  was  formulated  such  that  the  viscous 

*  I 

stresses  {o>vls  damp  out  the  highest  frequency  of  an  element  by 
approximately  the  specified  fraction  of  critical  element  damping  (v). 
The  corotational  viscous  normal  stresses  ^  were  calculated  in 

SAMS0N2  using  Equations  4.61,  4.62,  and  4.63. 


-11  ' 


=  {P¥13}  +  (s;p 


:4.6i) 


where. 


{PV1S} 

^!s> 


hydrostatic  viscous  normal  stresses,  and 
deviatoric  viscous  normal  stresses. 


{Pvis} 


2uB 

C<J. 


max 


/dll  +  d22  +  d33  s 

( - 3 - ) 


where, 

u 

B 


max 

AAA 

dH’d22,d33 


(4.62) 


specified  fraction  of  critical  element  damping, 

bulk  modulus  computed  in  the  stress  routines, 

maximum  frequency  of  the  element,  and 

normal  corotational  velocity  strains  in  the  x, 
y,  and  2  directions  that  were  calculated  in 
Section  4.2.3. 
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2uG 

“max 


4k  A 

+ 

22  33 
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(4.62) 


where, 

G  =  shear  modulus  for  a  given  material  model,  and 
d„  *  corotational  velocity  normal  strains,  1«1,  2,  3. 

*  y  j  $ 

The  corresponding  corotational  viscous  shear  stress  t  was  calculated 
in  SAMS0N2  using  Equation  4.64. 


;vis  „  hi  (4.64) 

“max  ^ 


The  development  for  Equations  4.61-4.64  and  the  equations  that  were  used 
in  order  to  approximate  u>  are  not  shown  here,  but  can  be  found  in 


Reference  1.  The  total  corotational  stresses  for  the  current  time  step, 
including  the  viscous  stresses  that  were  calculated  in  Equations  4.61- 
4.64,  were  calculated  in  SAMS0N2  using  the  expressions  shown  in 
Equation  4.65,  where  the  first  terms  on  the  righthand  side  of  the 
expressions  were  calculated  previously  in  the  corresponding  stress 
routines.  These  expressions  were  calculated  in  subroutine  SDAMP  and 
were  also  verified  through  comparisons  with  References  1  and  16. 
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These  modified  total  stresses  were  then  transformed  in  the  SAMS0N2  code 
from  the  corotational  coordinate  system  (x,y)  to  the  global  coordinate 
system  (x,y)  using  Equations  4.66  and  4.67. 


cos2© 

sin2© 

2sin  0  cos  © 


sin20  -sin  0  cos  0 

cos20  sin  0  cos  0 

-2sin  ©  cos  0  cos2©  -sin20 


(4.66) 


where. 


Tn  *  transformation  matrix  used  in  SAMS0N2  in  order  to  convert 

0 

the  corotational  stresses  (o)  to  stresses  in  the  global 
coordinate  system  (a),  and 

e  *  current  total  rotation  which  has  occurred  at  an  element 
quadrature  point  as  computed  by  Equation  4.13. 


{ c }  =  global  components  of  the  stresses  that  were  computed 

in  Equation  4.65,  and 

(o}n+1  =  stress  values  that  were  computed  using 

Equation  4.65. 

An  error  was  discovered  in  the  stress  transformation  matrix,  [T  1,  when 
comparisons  were  made  with  References  1,  7,  and  14.  The  correct 
transformation  matrix  which  should  be  used  in  the  SAMS0N2  formulation  is 
shown  in  Equation  4.68. 

cos2e  sin2e  -2sin  0  cos  9 

T  a  [R]T  =  $ln2e  cos2e  2sin  e  cos  e  (4.68) 

6 

sin  e  cos  9  -sin  e  cos  e  cos2e  -sin2e 

where, 

[R]T  *  transpose  of  the  strain  transformation  matrix  calculated 
in  Equation  4.14. 

It  should  be  mentioned  that  the  incorrect  transformation  matrix 
(Equation  4.66)  was  used  in  the  formulations  for  the  three-node 
triangular  (3NT),  4NQ,  5NT,  6NT,  and  8NQ  isoparametric  continuum  finite 
elements.  This  error  can  be  corrected  by  replacing  the  FORTRAN 
statement 

CALL  STPRO  (temp2, shat, stress) 

which  is  formulated  in  the  VnFRCN  subroutines  (n«3,4,5,6,  and  8)  with 
the  following  FORTRAN  statement: 

CALL  GTPRD  (temp2, shat, stress, 3, 3,1). 


The  internal  nodal  force  equations  which  were  formulated  into  the 
SAMS0N2  code  were  found  to  be  correct  when  compared  to  References  1,  9, 
13,  and  16.  However,  one  error  which  was  related  to  the  internal  force 
calculation  was  found.  This  error  was  in  the  formulation  for  the  stress 
transformation  matrix.  This  is  a  major  error  and  needs  to  be  corrected 
in  the  current  versions  of  SAMS0N2  code  in  order  for  the  program  to 
provide  correct  results. 

4. 2. 4. 6  Evaluation  of  the  Internal  Strain  Energy 

The  equations  that  were  formulated  in  the  SAMS0N2  code  for  the 
determination  of  the  internal  strain  energy  are  presented  in  this 
section.  Strain  energy  is  a  function  of  the  Cauchy  stresses  and  the 
velocity  strains  within  an  element.  The  internal  strain  energy  is 
calculated  correctly  for  these  conjugate  measures  of  strain  rate  and 
stress. 

The  basic  relation  which  is  used  for  determining  the  internal 
strain  energy  is  shown  in  Equation  4.69. 

U  *  i  f  {e}T  { o }  dV  (4.69) 

v 

where, 

U  *  total  internal  strain  energy  within  a  system,  and 

{e},{o}  =  stresses  and  strains  within  a  system. 


However,  this  equation  was  altered  in  order  to  be  consistent  with  the 
SAMS0N2  formulation.  The  resulting  expressions  which  were  used  in  the 
SAMS0N2  code  for  the  determination  of  the  total  internal  strain  energy 
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are  shown  in  Equations  4.70  and  4.71.  It  should  be  noted  that  a  2*  term 
was  not  included  in  the  determination  of  the  internal  strain  energy  for 
axisymmetric  analyses  because  it  is  canceled  in  the  determination  of  the 


energy  error. 


Un+1  =  un 


JE  n  n 


+  i  E  {  E  E  [rah  .h .  |J(Ci  ,n-  )l  Ae]l 
1=1  i=l  j=l  1  J  11 


(4.7C) 


where. 


un+1,un 


total  internal  strain  energies  of  the  system  for  the 
current  and  previous  time  steps,  respectively, 

total  number  of  elements  in  the  system,  and 

increment  of  internal  strain  energy  computed  for  an 
elemental  integration  point  as  shown  in 
Equation  4.71. 


nstres 


Ae  »  E  UeU^  ,nj)}n+*  x({a(c^  ,nj))n+1  +  {oU.j,nj)}n)  (4.71) 


where. 


nstres  =  3  for  plane  stress  analyses,  and 

=  4  for  plane  strain  or  axisymmetric  analyses; 


I 


4eU.-,n  )  *  incremental  strains  that  were  computed  in 

J  Equation  4.16  of  Section  4.2.3  for  a  particular 

Gaussian  quadrature  point;  and 

MSi,n.)}n  *  total  stresses  that  were  computed  using  the 
J  appropriate  material  law  routines  for  the  c 

and  previous  time  steps,  respectively. 


current 


Equations  4.70  and  4.71  were  verified  through  comparisons  with 
Reference  1.  The  value  for  the  total  internal  strain  energy,  Un+1,  was 
stored  in  the  strs  array  for  use  later  in  the  energy  error  solution. 
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4.2.5  Determination  of  the  External  Nodal  Forces 


and  the  Corresponding  External  Work 
The  calculations  that  were  involved  in  the  determination  of  the 
external  nodal  forces,  {Fgxt},  and  the  external  work  (W)  are  discussed 
in  this  section.  These  quantities  were  calculated  in  subroutines  FREEFD 
and  FREEF2  of  the  SAMS0N2  code.  Some  of  the  equations  that  were  used  in 
these  calculations  are  also  provided.  Formulation  errors  that  were 
discovered  in  these  two  subroutines  are  also  mentioned. 

The  SAMS0N2  code  contains  the  following  five  different  load  types 
(see  p.  50  in  Reference  3): 

1)  axisymmetric  pressure  load  line 

2)  plane  pressure  load  line 

3)  initial  impulse  load  line 

4)  force  line 

5)  displacement  history  load  line 

with  a  load  line  referring  to  the  consecutive  nodes  along  which  the  load 
was  applied.  The  external  nodal  forces  were  computed  in  accordance  with 
the  particular  load  type  that  was  used  in  the  problem  solution. 
Calculations  of  the  external  nodal  forces  were  done  in  approximately  the 
same  manner  for  each  of  the  five  load  types.  These  calculations  are 
discussed  for  each  of  the  load  types  in  the  following  paragraphs. 

Initial  nodal  velocities  were  calculated  for  the  initial  impulse 
load  line.  These  initial  nodal  velocities  were  calculated  prior  to  the 
start  of  the  SAMS0N2  explicit  integratioi.  *:heme.  These  initial  nodal 
velocities  were  calculated  in  subroutine  FREEFD  of  the  SAMS0N2  code 
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using  Equation  4.72  for  the  case  in  which  the  initial  impulse  load  was 


input  in  terms  of  the  normal  and  tangential  components. 
Vx  =  ^Vt  cos  0  "  vn  sin  0) 


V„  =  i(\L  sin  0  +  V„  cos  0) 


(4.72) 


where, 


VVy 


VfVn 


initial  nodal  velocities  transformed  from  the  input 
normal  and  tangential  components  to  the  global  x  and 
y  components, 

tangential  and  normal  components  for  the  input 
initial  impulse  load,  and 

angle  shown  in  Figure  4.4. 


The  sign  convention  that  was  used  with  this  impulse  load  is  illustrated 
in  Figure  4.4  with  the  normal  and  tangential  components  (Vn  and  Vt) 
shown  acting  in  the  positive  directions. 


Figure  4.4  Initial  Nodal  Velocities  for  an  Impulse  Load  Input  in 
Terms  of  Normal  and  Tangential  Components. 


Equation  4.72  was  found  to  be  -fncorrect.  The  correct  formula  for 

determining  the  initial  nodal  velocities  (V  and  V  )  is  shown  in 

x  y 

Equation  4.73. 


vx  =  *(vt  cos  0  +  Vn  sin  e) 
Vy  =  i(Vt  sin  e  -  Vn  cos  0) 


(4.73) 


In  comparing  Equation  4.72  to  Equation  4.73,  the  only  differences 

between  the  two  equations  are  the  signs  assigned  to  the  normal  velocity 

component  terms.  For  the  case  in  which  the  impulse  load  was  input  in 

global  x-  and  y-components ,  the  values  for  V  and  V  were  simply  set 

x  y 

equal  to  these  input  values.  The  values  for  V  and  V  were  stored  in 

x  y 

their  respective  locations  within  the  velocity  (V)  array  for  each 
segment  along  the  load  line.  Therefore,  all  of  the  nodes,  except  for 
the  first  and  last  nodes  on  the  load  line,  had  contributions  from  both 

of  the  segments  on  the  load  line  of  which  they  were  a  part.  Two  node 

factors,  factors  by  which  the  nodal  velocities  determined  in 

Equation  4.72  were  multiplied  (values  usually  equal  2.0),  were  used  in 
order  to  increase  the  values  for  the  Initial  velocities  for  these  first 
and  last  nodes  such  that  they  matched  the  values  for  the  other  nodes  on 
the  load  line.  Upon  determination  of  these  initial  nodal  velocities,  no 
more  calculations  were  executed  for  the  Initial  impulse  load  and 

subroutine  FREEFD  was  never  used  again. 

The  calculations  that  were  used  In  order  to  determine  the  external 
nodal  forces  (nodal  velocities  and  displacements)  for  the  displacement 
history  load  line  were  similar  to  those  for  the  impulse  load  line. 
However,  an  input  load  line  history  (displacement  versus  time  curve)  was 
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used.  Therefore,  the  first  calculation  that  was  performed  was  the 
determination  of  the  displacement  value  for  the  current  time  step  or  the 
total  solution  time.  This  displacement  value,  un+1,  was  calculated  in 
subroutine  PRESS.  This  displacement  value  was  then  stored  in  the 
locations  within  the  nodal  displacement  array  (xl)  corresponding  to  the 
specified  component  (x  or  y)  for  the  displacement  history  load  line. 
Nodal  velocities  were  also  computed  using  Equation  4.74  and  stored  in 
the  appropriate  locations  within  the  velocity  array  (V).  No  other 
calculations  were  performed  after  the  determination  of  these  external 


forces  (nodal  velocities  and  displacement)  in  the  FREEF2  subroutine. 


,n+l  _  un+1  -  un 


(4.74) 


where. 


n+1  n 
u  ,u  ■ 


nodal  velocity  computed  for  the  current  time  step 
£.t,  and 

displacement  values  obtained  from  the  load  line 
history  data  for  the  current  and  previous  time 
steps,  respectively. 


The  values  for  the  external  nodal  forces,  (Fext>,  that  were 
computed  for  the  force  line  load  were  obtained  using  the  subroutine 
PRESS  in  the  SAMS0N2  code.  The  values  for  the  external  nodal  forces 
were  set  equal  to  the  force  value,  Pn+1,  which  was  obtained  from  the 
force  line  history  data  for  the  current  time  step  using  PRESS.  These 
values  were  then  stored  in  the  fored  array  in  the  appropriate  locations 
corresponding  to  which  force  component  (x  or  y)  of  the  load  line  was 
specified  in  the  input  data.  Two  other  parameters  were  computed  for  the 
force  line  load.  These  parameters  were  the  total  external  work,  Wn+*, 


v  *- 


V 


and  the  total  linear  impulse,  WIM  ,  for  the  current  time  step.  These 


two  parameters  were  calculated  using  Equations  4.75  and  4.76. 

,  , ,  ndnod  •  , . 

Wn+1  =  VJn  +  J4t(Pn  1  +  Pn)(  r  u"  *) 


(4.75) 


where. 


w"*V 


p"*V 


ndnod 


total  external  work  computed  for  a  particular  force 
line  for  the  current  and  previous  time  steps, 

values  for  the  forces  for  the  current  and  previous 
time  steps  that  were  obtained  from  the  force  line 
history  data  using  PRESS, 

number  of  nodes  on  the  force  line,  and 


nodal  velocity  values  for  the  nodes  on  the  load 
line. 


WIMn+1  *  WIMn  +  *At(Pn+1  +  Pn}( ndnod) 


(4.76) 


where, 


WIMn+1 ,WIM 


*  total  linear  Impulse  for  the  externally  applied 
force  for  the  current  and  previous  time  steps. 


These  equations  were  verified  through  comparisons  with  References  1,  5, 
9,  and  16.  The  values  for  Wn+1,  Pn+1,  and  WIM1*-1  were  all  calculated  in 
FREEF2  and  were  then  stored  in  the  strs  array. 

The  calculations  for  the  external  nodal  forces  for  the  plane  and 
axi symmetric  pressure  load  lines  were  very  similar  to  those  for  the 
force  line  load.  However,  the  equations  that  were  used  in  the 
calculations  Included  the  areas  over  which  the  pressures  were  being 
applied.  The  expressions  that  were  used  in  order  to  determine  the 
external  nodal  forces  for  the  plane  and  axlsymmetric  pressure  load  lines 
are  shown  in  Equation  4.77  and  4.78. 
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f  =  iPn+1  t(yc2  -  ycl) 

n+l  (4.77) 

f  »  iP  1  t(xc2  -  xcl) 

where, 

f  ,f  2  external  nodal  force  components  that  were  computed 

y  for  a  particular  segment  along  a  plane  pressure  load 

line, 

Pn+*  2  current  pressure  value  that  was  obtained  from  the 

input  load  line  history  data  using  PRESS, 

t  2  thickness  of  the  elements,  and 

xcl, ycl  2  current  displaced  coordinates  for  node  I  of  the  load 

line. 


fx  »  *Pn+1  r(yc2  -  ycl) 
f  2  *Pn+1  r(xc2  -  xcl) 
where. 


(4.78) 


f  ,f  2  external  nodal  force  components  that  were  computed 
x  y  for  a  particular  segment  along  an  axisymmetric 

pressure  load  line,  and 


r  «  i(xci*xc2)  -  x-distance  from  the  axis  of  symmetry  to  the 

middle  of  the  load  line  segment. 


Figure  4.5  shows  an  example  segment  from  a  pressure  load  line  with  the 

applied  pressure  shown  acting  in  the  positive  direction.  The  external 

nodal  force  components  (f  and  f  )  were  stored  in  the  appropriate 

x  y 

locations  within  the  forcd  (external  nodal  force)  array  that 
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Figure  4.5  Typical  Segment  Along  a  Pressure  Load  Line. 


Equations  4.79  and  4.80  are  the  expressions  that  were  used  in  subroutine 
FREEF2  of  the  SAMS0N2  code  in  order  to  determine  the  current  total 
external  work  for  the  plane  and  ax i symmetric  pressure  load  lines, 
respectively. 


Wn+1  -  Wn  ♦  i(Pn+1  ♦  P)[(pvol"+1  -  pvoln)] 


(4.79) 


where , 
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-  tuyj ) (uxj^j)0 


(4.80) 


xcj.ycj  * 


parameter  calculated  for  time  step  n  which 
incorporates  the  total  nodal  displacement  components 
and  the  x-  and  y-components  of  the  distances  between 
two  adjacent  nodes  on  the  load  line  (area  units), 

the  originally  input  x-  and  y-coordinates  for 
node  I , 
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u  T,u  T  =  total  displacement  components  for  node  I  in  the 

X1  yl  mesh, 

ndn  =  ndnod  -  1  =  number  of  nodes  on  a  load  line  -  1,  and 

r°  =  t  for  plane  analyses 

=  $(xcl  +  xc2)  for  axisynmetric  analyses. 

However,  in  a  preliminary  comparison  made  between  these  equations  and 
Equation  4.81  (obtained  from  References  1,  9,  and  16),  the  external  work 
formulation  in  the  SAMS0N2  code  for  the  pressure  line  loads  was  not 
verified. 

Wn+1  =  Wn  +  i4t{un+i>T((F2xt>  +  (f£*J))  (4.81) 

The  linear  impulse  for  the  external  nodal  forces  was  obtained  using 
Equation  4.82.  This  equation  was  found  to  be  correct  when  compared  to 
material  in  Reference  5. 

_4.i  n  ndn  n+1  ndn  .  .  pn+l  pn 

WIMn  1  =  WIMn  +  i&t[(  Z  fj+l)  +  (  Z  f"+1)]*  —  — t—  (4.82) 

1-1  x  1-1  y  Pn  1 

The  pvoln+1,  Pn+1,  Wn+1,  and  WIMn+*  parameters  were  all  stored  in  the 
strs  array  for  use  in  the  calculations  for  the  next  time  step. 

The  majority  of  the  formulation  in  subroutines  FREEF2  and  FREEFD 
that  was  used  in  the  determination  of  the  external  nodal  forces, 
external  work,  and  the  linear  impulse  was  found  to  be  correct.  However, 
two  errors  in  the  formulation  were  found  and  need  to  be  corrected  in 
order  to  obtain  accurate  results  for  all  analyses.  These  errors  were  in 
the  determination  oi  the  initial  nodal  velocities  for  the  impulse  load 
and  the  calculation  of  the  external  work  ^or  the  pressure  load  lines. 
One  additional  note  regarding  the  equations  for  the  axisymmetric 
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pressure  load  is  that  the  factor  2*  was  not  included  in  any  of  the 
equations.  This  elimination  of  the  factor  2*  was  consistent  with 
previous  axisymmetric  formulations. 


4.2.6  Solution  to  the  Equation  of  Motion 

This  section  contains  a  brief  discussion  regarding  the  general 
equation  of  motion  followed  by  a  discussion  to  the  solution  of  the 
equation  of  motion  formulated  in  the  SAMS0N2  code.  The  equation  of 
motion  is  a  mathematical  expression  through  which  the  dynamic 
displacement  values  for  acceleration,  velocity,  and  displacement  are 
defined.  The  general  expression  for  the  equation  of  motion  is  as 


follows: 


[M](u>  +  [C](u>  +  [K]{u)  -  P(t) 


(4.83) 


where. 


(uMul.tu) 


mass  matrix  for  the  system, 

viscous  damping  matrix  for  the  system, 

system  stiffness  matrix, 

dynamic  displacement  values  for  the 
acceleration,  velocity,  and  displacement 
occurring  within  a  system,  and 

time  varying  externally  applied  load. 


Each  term  in  Equation  4.83  corresponds  to  a  type  of  force  with  the  terms 
on  the  lefthand  side  of  the  equation  being  the  Inertial,  damping,  and 
elastic  forces,  respectively.  The  solution  to  the  equation  of  motion 
provides  the  values  for  the  accelerations,  velocities,  and  displacements 
for  a  particular  point  in  time.  And,  when  the  solution  to  the  equation 


re- 


of  motion  is  obtained  for  different  points  in  time,  the  time  histories 
for  these  three  parameters  can  be  obtained. 

The  equation  of  motion  that  was  formulated  in  subroutine  SOLVE  of 
the  SAMS0N2  code  for  time  step  n  is  shown  in  Equation  4.84. 

[M]('u)n  *  [C]^"  *  <Ffnt)n  =  (Fext)n  (4.84) 

where, 

[M]  =  lumped  diagonal  mass  matrix  that  was  determined  in 

Section  4.2.2, 

[C]  =  damping  matrix  shown  in  Equation  4.85, 

(F.  )  =  internal  nodal  forces  that  were  computed  in 

Section  4. 2. 4. 5,  and 

(F  }  =  external  nodal  forces  that  were  calculated  in 

ext  Section  4.2.5. 

This  equation  required  the  following  two  initial  conditions:  (1)  { u>°  * 
(u(0)}  and  (2)  (u)’*  ■  (u(O)h  The  accuracy  of  the  solution  due  to 
these  initial  velocity  conditions  is  affected  only  slightly.  The 
damping  matrix  that  was  used  in  the  SAMS0N2  formulation  is  of  the  form 
shown  in  Equation  4.85. 

[C]  •  Cj CM3  ♦  C2[K]  (4.85) 

where, 

Cj  ■  mass  proportional  damping  factor,  and 
C2  s  stiffness  proportional  damping  factor. 

The  stiffness  proportional  damping  term,  was  included  in  the 

determination  of  the  nodal  internal  forces  through  the  use  of  linear 
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artificial  viscosity.  The  mass  proportional  damping  term  was 
substituted  into  Equation  4.84  for  [C]  which  resulted  in  Equation  4.86. 


[M]{u'}n  ♦  CjCMKul"  *  {fint}n  *  /Fextln 


(4.86) 


The  mass  proportional  damping  term  is  only  used  for  dynamic  relaxation 
problems  (see  Reference  1)  in  order  to  simulate  static  equilibrium 
through  the  use  of  sufficient  damping.  The  value  for  is  determined 
with  the  use  of  an  estimate  of  the  minimum  frequency  of  a  system  as 
shown  on  pp.  128-132  of  Reference  3.  The  values  for  the  nodal 
accelerations  at  ar>  time  step  were  obtained  from  the  solution  to 
Equation  4.86.  Equation  4.8?  shows  the  expression  that  was  used  in 
SAMSCN2  in  order  to  determine  the  values  for  the  nodal  accelerations  for 
time  step  n,  <uin. 


(4.87) 


where. 


[M]“  «  inverse  of  the  diagonal  lumped  mass  matrix,  and 

{u}n~*  «  nodal  velocities  at  the  half  time  step. 

Equation  4.87  was  obtained  by  substituting  the  expression  shown  in 
Equation  4.88  into  Equation  4.86  and  rearranging  the  terms. 


f»,n  at  ", n  A  ,*,n-J 
{u}  »  (u)  +  {u} 


(4.88) 


Equation  4.88  was  obtained  through  the  use  of  central  finite  difference 
expressions.  The  inversion  of  the  diagonal  lumped  mass  matrix,  [M] 
in  Equation  4.87  is  trivial  and  only  involves  the  Inversion  of  each  term 
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on  the  diagonal.  It  should  be  noted  that  the  acceleration  values  that 
are  output  for  a  particular  analysis  involving  dynamic  relaxation  are 
not  equal  to  the  values  obtained  from  Equation  4.87.  Instead,  the 
values  that  are  output  for  the  accelerations  are  obtained  from  the 
solution  to  Equation  4.89.  These  values  are  the  actual  accelerations  in 
a  dynamic  analysis. 

(u}n  =  CM]'1({Fext}n  -  <Fint>n)  (4.89) 

The  values  for  the  nodal  displacements  and  velocities  were  calculated  in 
SAMS0N2  using  the  central  difference  expressions  shown  in  Equations  4.90 
and  4.91. 

(G}n+*  *  (u}n‘*  +  At  (u)n  (4.90) 

where, 

{ u }  4  3  nodal  velocities  for  the  previous  time  step,  and 

(u}n  *  values  for  the  nodal  accelerations  that  were 

calculated  using  Equation  4.87. 

(U)n+1  «  (U)n  ♦  At(u)n+i  (4.91) 

where, 

(u)n+1,{u}n  *  nodal  displacements  for  the  current  and  previous 
time  steps. 

The  values  for  the  accelerations,  velocities,  and  displacements  that 
were  calculated  using  Equations  4.87,  4.89,  4.90,  and  4.91  were  adjusted 
in  the  SOLVE  routine  to  reflect  the  boundary  conditions  of  the  problem. 
All  of  the  equations  in  the  SOLVE  routine  were  verified  through 
comparison  with  References  1,  9,  and  16.  It  should  be  noted  that  the  ?* 


factor  does  cancel  out  in  the  equation  of  motion  for  axisymmetric 
analyses. 

The  values  (u,  u,  and  u)  that  were  obtained  by  solving 
Equations  4.87,  4.89,  4.90  and  4.91  were  stored  in  the  nodal 

acceleration  (a),  velocity  (v),  and  displacement  (xl)  arrays, 
respectively,  for  use  in  the  next  integration  cycle  as  shown  in 
Table  4.8.  This  table  shows  the  order  for  which  the  calculations  that 
were  exhibited  in  Section  4.2  were  performed  in  the  SAMS0N2  explicit 
integration  scheme. 


4.3  Summary  of  the  SAMS0N2  Formulation  Errors 


The  8NQ  higher-order  finite  element  formulation  in  the  SAMS0N2  code 
was  investigated  and  verified.  Formulation  errors  were  discovered  in 
SAM SON 2  through  this  Investigation  and  verification  process.  The 
following  is  a  list  of  errors  and  possible  errors  that  were  discovered 
in  the  SAMS0N2  formulation: 

(1)  Twelve  Gauss-legendre  abscissae  and  weight  coefficients  were 
incorrectly  typed  Into  the  GAUSS1  subroutine  (see  Table  4.3). 

(2)  The  equation  that  was  used  in  subroutine  V8ASME  in  order  to 
determine  the  mass  for  an  axisymmetric  8NQ  isoparametric  continuum 
element  is  incorrect  (see  Equation  4.7  and  the  associated 
discussion).  Also,  a  possible  approximation  error  exists  related 
to  the  use  of  Equation  4.7  in  the  determination  of  mass  for  both 
the  4NQ  and  8NQ  elements. 
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(3)  A  possible  error  exists  in  the  formulation  for  the  strain 
transformation  matrix  (Equation  4.14).  The  term  sin  0  is  set  equal 
to  zero  in  subroutine  VSTRAN  for  |  0  |  <  1.0  x  10“^. 

(4)  An  error  may  result  if  the  biaxial  elastic-plastic  plane  stress  and 
biaxial  elastic-plastic  plane  strain  or  axisynmetric  material  laws 
are  used  in  analyses  involving  improper  materials.  The  von  Mises 
failure  criterion  is  only  valid  for  ductile  or  semiductile  metals. 

(5)  The  equation  in  subroutine  STRES3  that  is  used  in  order  to 
determine  the  strain  increment  in  the  z-direction  (Ae^)  is 
incorrect  for  plane  stress  analyses  where  yielding  has  occurred 
(see  Equations  4.32  and  4.33  and  corresponding  discussion). 

(6)  The  formulation  for  the  AFWL  "engineering"  model  in  subroutine 
STRES9  (see  Section  4. 2. 4. 4)  contains  the  following  four  errors: 

(a)  The  out-of-plane  stress  (oz)  is  set  equal  to  zero  for  plane 
strain  analyses  which  contain  elements  having  a  thickness  not 
equal  to  1.0. 

(b)  The  factor  Fa  which  is  the  portion  of  the  strain  which  reloads 
out  of  the  tension  cut-off  was  not  Included  in  the  STRES9 
formulation. 

(c)  A  possible  error  exists  if  the  computed  value  for  the 
volumetric  strain  (evo1)  exceeds  the  last  limiting  volumetric 
strain  value  on  the  hydrostat  curve.  This  error  would  occur 
if  the  bulk  modulus  for  the  last  loading  segment  on  the 
hydrostat  curve  was  not  equal  to  the  lockup  bulk  modulus. 

(d)  Another  possible  error  may  occur  if  the  value  that  is 
calculated  for  the  current  hydrostatic  pressure  (Pn+*)  exceeds 
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the  value  for  the  hydrostatic  pressure  which  is  associated 
with  the  last  input  yield  function  point.  If  this  condition 
occurs,  then  the  value  for  YJ  is  set  equal  to  the  last  yield 
point  value  which  may  or  may  not  be  the  correct  value. 

(7)  The  transformation  matrix  which  is  used  in  order  to  transform  the 
stress  components  from  the  corotational  coordinate  system  to  the 
global  coordinate  system  is  incorrect  (see  Equations  4.66  and  4.68 
and  associated  discussion).  This  error  exists  for  all  of  the  plane 
and  axisymmetric  continuum  elements. 

(8)  The  equations  in  subroutine  FREEFD  that  are  used  to  calculate  the 
initial  nodal  velocities  for  an  initial  impulse  load  which  has  been 
input  in  normal  and  tangential  components  are  incorrect  (see 
Equations  4.72  and  corresponding  discussion). 

(9)  The  external  work  equations  in  subroutine  FREEF2  that  are  used  in 
order  to  compute  the  external  work  due  to  a  plane  or  axisymmetric 
pressure  load  line  appear  to  be  incorrect  based  on  a  preliminary 
investigation  that  was  performed  (see  Equations  4.79,  4.80,  and 
4.81  and  corresponding  discussion). 

(10)  The  column  labels  that  are  associated  with  the  output  from  the  AFWl 
"engineering"  model  are  incorrect.  The  current  hydrostatic 
pressure  and  the  current  volumetric  strain  are  labeled  sig-yld  and 
rot,  respectively. 

These  errors  in  the  current  finite  element  formulation  for  the  SAMS0N2 

code  need  to  be  corrected  an.1  tested. 


CHAPTER  5 


Effect  of  the  Corrected  Formulation  on  the  Finite  Element  Analyses 

This  chapter  contains  a  discussion  of  the  changes  that  occurred  in 
the  solutions  to  the  problems  from  Chapter  3  after  the  corrections  were 
made  to  the  SAMS0N2  finite  element  formulation  based  on  the  results  from 
Chapter  4.  The  purpose  of  this  chapter  was  to  determine  if  more 
consistent  and  more  reliable  results  could  be  obtained  through  the  use 
of  the  corrected  finite  element  formulation.  It  should  be  noted  that 
some  of  the  errors  that  were  discovered  in  the  SAMS0N2  finite  element 
formulation  were  unrelated  to  the  problems  from  Chapter  3.  Hence,  the 
problem  solutions  were  unaffected  and  the  significance  of  these  errors 
could  not  be  determined. 

The  majority  of  the  corrections  that  were  made  to  the  SAMS0N2 
formulation  had  little  or  no  effect  on  the  results  for  the  problems  in 
Chapter  3.  For  Instance,  the  corrections  for  the  Gauss-Legendre 
abscissae  and  weight  coefficients  In  subroutine  GAUSS1  resulted  In  no 
change  for  any  of  the  8NQ  element  solutions.  However,  all  of  these 
solutions  used  an  order  of  integration  of  2.0.  These  corrected 
coefficient  values  may  be  significant  for  problems  in  which  a  higher 
order  of  integration  is  specified.  The  approximation  for  sin  e  that  is 
used  in  subroutine  VSTRAN  for  J  0  |  <  1.0  x  10'^  was  changed  to  sin 
9  *  sin  0  for  all  values  of  ©.  This  correction  was  very  insignificant 
as  might  be  expected.  The  equation  in  subroutine  STRES3  that  is  used  to 
calculate  the  incremental  strain  in  the  z-direction  (4e2)  was  corrected 
using  Equation  4.33.  This  correction  resulted  in  only  a  1%  difference 
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in  the  values  for  the  total  strains  (e  )  between  the  two  solutions.  No 
other  values  (displacements,  strains,  or  stresses)  were  affected  by  this 
correction  to  the  STRES3  routine.  The  correction  to  the  stress 
transformation  matrix  also  was  insignificant  in  its  effect  on  the 
problem  solutions.  Only  slight  changes  occurred  in  the  solutions  due  to 
this  modification.  The  use  of  the  r  approximation  in  the  determination 
of  the  diagonal  lumped  mass  matrix  for  an  axisymmetric  continuum  was 
found  to  have  no  effect  on  the  results  for  the  axisymmetric  problems. 
However,  the  equation  for  determining  r  was  in  error  in  the  original 
SAMS0N2  formulation.  The  correct  expression  for  r  (Equation  4.7)  was 
substituted  into  the  updated  formulation  which  resulted  in  significant 
changes  in  the  8NQ  axisymmetric  solutions.  This  corrected  r  expression 
resulted  in  very  significant  changes  in  the  8NQ  element  results  for  both 
the  ax i symmetric  wave  propagation  analysis  and  the  axlsymnetric  analysis 
of  the  soil -structure  interaction  problem,  figures  5.1,  5.2,  and  5.3 
show  the  relative  improvements  In  the  results  for  the  axisyemetrlc  wave 
propagation  when  compared  to  Figures  3.23,  3.24,  and  3.25.  As  shown, 
the  8NQ  solution  which  was  obtained  using  the  updated  formulation  no 
longer  lags  behind  the  4NQ  element  solution.  The  corrected  r 
formulation  also  improved  the  results  for  the  axisymmetric  analysis  of 
the  soil -structure  interaction  problem  which  was  obtained  using  the  8NQ 
isoparametric  element.  Table  5.1  when  compared  to  Table  3.6  shows  the 
relative  improvement  between  the  two  8NQ  solutions.  The  8NQ  solution 
which  was  obtained  using  the  corrected  r  formulation  exhibited  a  much 
better  agreement  with  the  4NQ  solution.  A  significant  improvement  in 
the  8NQ  solution  to  the  soil-structure  interaction  problem  for  the 
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Table  5.1  Comparison  of  Y-Displacement  Values  Between  the  4NQ  and 
3NQ  Solutions  for  the  8NQ-6  Element  Structural  Mesh  and 
the  4NQ-12  Element  Structural  Mesh  (1.0  psi  Peak  Load) 
(Soil-Structure  Interaction  Problem)  After  Mass 
Formulation  Corrected. 


Original 

Node  Location 

x-coordinate  y-coordinate 
(in)  (in) 

8NQ 

Y-Oisplacement 
t=.8xl0"3  sec 
(in  x  10'4) 

4NQ 

Y-Displacement 
t=.8xl0"3  sec 
(in  x  10’4) 

Percent 

Difference 

Based  on 

4NQ  Solution 

0.0 

72.0 

-.2231 

-.2112 

5.6% 

3.0 

72.0 

-.2068 

— 

— 

6.0 

72.0 

-.1961 

-.2009 

2.4% 

9.0 

72.0 

-.1882 

— 

— 

12.0 

72.0 

-.2154 

-.1972 

9.2% 

0.0 

66.0 

-.2076 

-.2073 

0.1% 

6.0 

66.0 

-.1894 

-.1961 

3.4% 

12.0 

66.0 

-.1971 

-.1912 

3.1% 

0.0 

60.0 

-.1890 

-.1996 

5.3% 

3.0 

60.0 

-.1869 

—  —  • 

— 

6.0 

60.0 

-.1916 

-.1878 

2.0% 

9.0 

60.0 

-.1941 

— 

•  •  • 

12.0 

60.0 

-.1827 

-.1831 

0.2% 

0.0 

54.0 

-.1886 

-.1949 

3.2% 

6.0 

54.0 

-.1882 

-.1785 

5.4% 

12.0 

54.0 

-.1777 

-.1710 

3.9% 

0.0 

48.0 

-.1731 

-.1898 

8.8% 

3.0 

48.0 

-.1795 

- — 

— 

6.0 

48.0 

-.1792 

-.1711 

4.7% 

9.0 

48.0 

-.1751 

— 

— 

12.0 

48.0 

-.1745 

-.1612 

8.3% 

0.0 

42.0 

-.1584 

-.1788 

11.41 

6.0 

42.0 

-.1685 

-.1639 

2.8% 

12.0 

42.0 

-.1724 

-.1563 

10.3% 
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applied  blast  load  was  also  observed.  The  remainder  of  the  errors, 
listed  in  Chapter  4  and  not  discussed  in  this  section,  were  not  related 
to  the  problem  solutions  in  Chapter  3. 

In  conclusion,  the  majority  of  the  SAMS0N2  formulation  errors  that 
were  presented  in  Chapter  4  were  insignificant  in  their  effect  on  the 
problem  solutions.  However,  the  corrected  r  formulation  rendered  a 
significant  improvement  in  the  8NQ  axisymmetric  analyses.  On  the  basis 
of  this  result,  the  corrected  SAMS0N2  formulation  does  provide  more 
consistent  and  more  reliable  results. 
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CHAPTER  6 


Conclusions  and  Recommendations 


The  primary  conclusion,  based  on  the  results  that  were  presented  in 
the  previous  chapters,  is  that  the  current  finite  element  formulation 
for  the  eight-node  quadrilateral  isoparametric  continuum  element 
provides  consistent  and  reliable  results  for  problems  which  involve 
plane  analyses.  However,  some  minor  errors  exist  in  the  general  SAMS0N2 
finite  element  formulation  that  need  to  be  corrected.  Also,  the  8NQ 
axisymmetric  continuum  formulation  contains  a  significant  error  in  the 
axisymmetric  mass  formulation  that  needs  to  be  corrected  in  the  current 
formulation  in  order  for  accurate  results  to  be  obtained.  Hence,  the 
errors  that  were  discovered  in  the  SAMS0N2  formulation  (see  Chapter  4) 
must  be  corrected  in  order  for  consistent  and  reliable  results  to  be 
obtained  for  both  8NQ  axisymmetric  analyses  and  general  continuum  finite 
element  analyses. 

Some  further  conclusions  that  were  drawn  based  on  the  results  from 
the  previous  chapters  are  as  follows: 

(1)  The  higher-order  element  (5NT,  6NT,  and  8HQ)  solution  schemes  are 
less  efficient  than  the  4NQ  solution  scheme.  This  conclusion  is 
based  on  the  fact  that  much  longer  execution  times  were  required 
for  the  higher-order  element  solutions. 

(2)  The  finite  element  formulation  for  the  6NT  Isoparametric  continuum 
element  works  properly  for  elastic  plane  analyses  of  flexure 
problems.  This  conclusion  is  based  on  the  good  correlation  in 
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results  that  was  obtained  in  the  cantilever  beam  analyses  between 
the  SAMS0N2  and  the  corresponding  analytical  solutions. 

(3)  The  5NT  isoparametric  continuum  elemert  should  only  be  used  as  a 
transition  element  for  which  it  was  designed.  This  conclusion  is 
based  on  the  5NT  results  for  the  cantilever  beam  problem. 

(4)  Smaller  time  step  values  are  necessary  for  solutions  which  involve 
higher-order  elements  as  opposed  to  the  values  that  are  required 
for  4NQ  element  solutions. 

(5)  The  higher-order  element  (6NT  and  8NQ)  solutions  were  more  accurate 
in  some  cases  compared  to  corresponding  4NQ  solutions  and  were 
obtained  using  fewer  nodes  and  elements. 

Reconmendations  for  future  work  are  presented  below.  These 

reconmendations  are  based  upon  both  the  work  that  has  already  been 

performed  and  the  results  that  were  obtained  during  this  study. 

(1)  The  8NQ  element  needs  to  be  used  with  the  corrected  formulation  in 
a  soil -structure  Interaction  (SSI)  problem  for  which  test  data  are 
available  in  order  to  determine  whether  the  8NQ  element  solution  is 
accurate  for  an  SSI  problem. 

(2)  Additional  problems  need  to  be  analyzed  using  the  8NQ  higher-order 
element  In  order  to  test  the  performance  of  the  SAMS0N2  8NQ  finite 
element  formulation  more  thoroughly. 

(3)  The  5NT  and  6NT  continuum  elements  need  to  be  investigated  more 
extensively  through  the  solutions  to  more  problems  and  the 
formulation  verifications. 

(4)  Problems  need  to  be  analyzed  using  discretizations  which  contain 
both  the  higher-order  and  the  4NQ  elements.  Results  for  these 
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problems  could  be  used  to  determine  whether  transitions  can  be  made 
between  the  different  elements  with  accurate  results  still  being 
obtained. 

(5)  The  corrected  finite  element  formulation  needs  to  be  further  tested 
since  some  of  the  errors  were  unrelated  to  the  problem  analyses 
that  were  performed  in  this  investigation.  For  instance,  the 
corrections  in  the  AFWL  "engineering"  model  need  to  be  further 
tested. 

(6)  A  further  investigation  and  verification  needs  to  be  undertaken 
regarding  the  equations  that  are  used  in  the  SAMS0N2  code  in  order 
to  compute  the  external  work  for  plane  and  axi symmetric  pressure 
loadings. 

(7)  If  the  8NQ  finite  element  formulation  produces  inaccurate  results 
in  future  problem  analyses,  then  the  following  two  items  should  be 
performed: 

(a)  development  of  alternative  schemes  for  determining  the  nodal 
masses,  and 

(b)  investigation  of  the  addresses  and  Indexes  that  are  used  in 
the  8NQ  finite  element  formulation  in  order  to  store  and 
recall  information  from  the  various  parameter  arrays. 


REFERENCES 


Belytschko,  T.  and  Robinson,  R.  R.,  SAMS0N2:  A  Nonlinear 
Two-Dimensional  Structure-Media  Interaction  Computer  Cocfe . 


r  Force  Weapons  Laboratory,  Kirtland  Air  Force 
Base,  New  Mexico,  January  1982. 


Schreyer,  H.  L.,  Personal  coronuni cation  with  regard  to  the 
higher-order  finite  element  formulation  in  the  SAMS0N2  code. 


Schreyer,  H.  1.,  et  al.,  SAMS0N2,  A  Nonlinear  Two-Dimensional 
Structure-Media  Interaction 


,  Air  Force  Weapons  Laboratory,  Kirtland  Air  Force 
Base,  New  Mexico,  June  1984. 


4.  Seemann,  D.  R.,  Personal  communication  with  regard  to  the 
higher-order  finite  element  formulation  in  the  SAMS0N2  code. 

5.  Clouqh,  R.  W.  and  Penzien,  J.,  Dynamics  of  Structures, 
McGraw-Hill,  New  York,  1975. 

6.  Berglund,  J.  W.  and  Rudeen,  D.  K.,  STEALTH  and  SAMS0N2 
Verification  Calculations,  NTE-TN-83-20,  Air  Force  Weapons 
Laboratory,  Kirtland  Air  Force  Base,  New  Mexico. 

7.  Higdon,  A.,  et  al . ,  Mechanics  of  Materials,  4th  ed.,  John  Wiley  & 
Sons,  New  York,  1985. 

8.  Paz,  M.,  Structural  Oynamics:  Theory  end  Computation,  2nd  ed..  Van 
Nostrand  Reinhold  Company,  New  York,  1985. 

9.  Belytschko,  T.  and  Hughes,  T.  J.  R.,  eds..  Computational  Methods 
for  Transient  Analysis,  Elsevier  Scienci  Publishers  B.V. , 
Amsterdam,  1$B3. 

10.  Fenves,  S.  J.,  et  al . ,  eds..  Numerical  and  Computer  Methods  in 
Structural  Mechanics,  Academic  Press,  New  York,  19/3. 

11.  Boresi,  A.  P.,  et  al..  Advanced  Mechanics  of  Materials,  3rd  ed., 
John  Wiley  &  Sons,  New  York,  1978. 

12.  Bathe,  K.  J.,  Finite  Slqynt  Procedures  in  Engineering  Analysis, 
Prentice-Hall,  Englewood  Cliffs,  New  Jersey,  1982. 

13.  Gallagher,  R.  H.,  Finite  Element  Analysis:  Fundamentals, 

Prentice-Hall,  Englewood  Cliffs,  New  Jersey,  1975. 

14.  Weaver,  W.  and  Johnston,  P.  R.,  Finite  Elements  for  Structural 
Analysis,  Prentice-Hall,  Englewood  Cliffs,  New  Jersey,  1984. 

15.  Zienklewicz ,  0.  C.,  The  Finite  Element  Method,  3rd  ed., 

McGraw-Hill,  New  York,  197$.  ~ 


16.  Belytschko,  T.,  Chiapetta,  R.  L.,  and  Bartel,  H.  D. ,  "Efficient 
Large  Scale  Non-Linear  Transient  Analysis  by  Finite  Elements." 
International  Journal  for  Numerical  Methods  in  Enqineerin 


1 


ELASTIC  CONT.  16  QUAD  ELTS  WITH  DISPLACEMENT/ FREE  RT  END 


27  16 

1 

6 

130  0.25 

1 

0 

3 

1 

980. 

CM  DYNE  DYNE 

1  6 

4 

ELASTIC  PLANE  STRESS 

0.0 

1.0 

1.0 

1.0 

0.0  0.0 

1 

9 

16.0 

10 

2.0 

18 

16.0 

2.0 

19 

0.0 

4.0 

27 

16.0 

4.0 

1  1 

2 

11 

10 

1 

8  8 

9 

18 

17 

1 

9  10 

11 

20 

19 

1 

16  17 

18 

27 

26 

1 

121 

1021 

1921 

901 

1801 

2701 

1  8  3 

1104X-DISPL  NODE-  1 
3104X-DZ5PL  NODE-  3 
5104X-DISPL  NODE-  5 
7104X-DISPL  NODE-  7 
14104X-DISPL  NODE-14 
3U4X-VEL  NODE-3 
5114X-VEL  NODE -5 
7114X-VEL  NODE-7 
2  S4 STRESS  AT  2 
4  S4STEZSS  AT  4 
«  V45TKSS  AT  6 
25  3 

SO  3 

2  15  3 

1  3  5  7  14 

2  4  4 

13  2  15 

0.  0.  2. 

8.  .01  10. 

14.  .02  35. 

1  15 


2 


10 

S.0E-5 

1.3831*2 

.02 


4.  2.53E-3 

12.  1.7071-2 


4.  6.17E-3 

14.  1.524E-2 


Figure  A.l 


One-Dimensional  Wave  Propagation  Problem:  Input  for  the 
4 NO  Element  Discretization. 


ELASTIC  CONT. 
69  16 

3 

980. 

1  6 


16  8-NODED  QUAD  ELTS  WITH  DISPLACEMENT/ FREE  RT  END 


CM  DYNE  DYNE 
ELASTIC  PLANE  STRESS 


1 

1 

9 

11 

3 

6 

10 

7 

2 

8 

S7 

65 

67 

59 

62 

66 

63 

58 

9 

3 

11 

13 

5 

7 

12 

8 

4 

16 

59 

121 

67 

69 

61 

1 

63 

68 

64 

60 

1  8 
1  8 
1  8 
1  8 


Figure  A. 2  One-Dimensional  Wave  Propagation  Problem:  Input  for  the 
8NQ  Element  Discretization. 


WSWW55 


W5 


W'V'V'V 


S 

tS 

rW1 


153 


1104X-DISPL  NODE *  1 
17104X-DISPL  NODE *17 
33104X-DISPL  NODE*33 
49104X-DISPL  N0DE*49 
35104X-DISPL  NODE*35 
17114X-VEL  NODE*17 
33114X-VEL  MODE* 3 3 
49114X-VEL  NODE-49 
2  S4STXZS5  El .2  *1 
2144STRESS  EL 2  #2 
2234SIRESS  EL2  #3 
2324STRESS  EL2  #4 
4  S4STRESS  EL4  #1 
4144STRESS  EL4  #2 
4234STRESS  EL4  *3 
4324STRESS  EL4  #4 
6  54STRESS  EL6  #1 
6144STRESS  EL6  92 
6234STRESS  EL6  #3 
63245TRESS  EL6  «4 


25 

3 

50 

3 

1 

5 

3 

17 

33 

49 

35 

4 

6 

5 

2 

1 

1 

10 

0. 

0. 

2, 

8.0E-5 

8. 

.01 

10. 

1.383E-2 

16. 

.02 

35. 

.02 

1  5 


4.  2.93E-3 

12.  1.707E-2 


6.  6.17E-3 

14.  1.924E-2 


Figure  A.  2 


One-Dimensional  Wave  Propagation  Problem:  Input  for  the 
8NQ  Element  Discretization  (Continued). 


STATIC  ANAL,.  OF  A  CANT.  BEAM  USING  64  4NQ  ELEMS .  6/20/8S  BY  S.  MILLER 


85 

64 

1 

5 

4000 

3.4E-5 

1 

0 

2 

0 

0 

0 

1 

32.2 

FT 

1 

3 

4 

4 

-NODE 

PLANE 

CONTINUUM 

ELEMENTS 

0.01 

2.0 

16.0 

4176 

.  0E6 

0.0 

0.00 

0.0 

1 

0.0 

0.0 

17 

16.0 

0.0 

18 

0.0 

1.0 

34 

16.0 

1.0 

35 

0.0 

2.0 

51 

16.0 

2.0 

52 

0.0 

3.0 

68 

16.0 

3.0 

69 

0.0 

4.0 

85 

16.0 

4.0 

1 

1 

2 

19 

18 

1 

16 

16 

17 

34 

33 

1 

17 

18 

19 

36 

35 

1 

32 

33 

34 

51 

50 

1 

33 

35 

36 

53 

52 

1 

48 

50 

51 

68 

67 

1 

49 

52 

53 

70 

69 

1 

64 

67 

68 

85 

84 

1 

111 


1811 

3511 

5211 

6911 

1  1  2 
51200 
1  50 
8  50 
5000  5000 

1  17  -1  13 

0.0  3.1E-3  100.0  1.0  100.0 

69 


Figure  A. 3 


Cantilever  Beam  Problem:  Input  for  the  4NQ  Element 
Discretization. 


J*  .* 


,TIC 

ANAL, 

,  OF 

A  CANT.  BEAM  USING  32  5 NT 

ELEMS . 

6/20/85  BY 

61 

32 

1 

5 

4000 

3.4E-5 

1 

0 

2 

1 

32.2 

FT 

1 

3 

4 

5 

-NODEC 

i  TRIANGULAR  CONTINUUM 

ELEMENTS 

2.0 

2.0 

16.0 

4176 

.  0E6 

0.0  0.0 

0, 

.0 

1 

0.0 

0.0 

9 

16.0 

0.0 

10 

0.0 

1.0 

26 

16.0 

1.0 

27 

0.0 

2.0 

35 

16.0 

2.0 

36 

0.0 

3.0 

52 

16.0 

3.0 

53 

0.0 

4.0 

61 

16.0 

4.0 

1 

1 

2 

27 

11 

10 

1 

2 

2 

3 

29 

14 

13 

1 

3 

3 

4 

29 

15 

14 

1 

4 

4 

5 

31 

18 

17 

1 

5 

5 

6 

31 

19 

18 

1 

6 

6 

7 

33 

22 

21 

1 

7 

7 

8 

33 

23 

22 

1 

8 

8 

9 

35 

26 

25 

1 

9 

28 

27 

2 

11 

12 

1 

10 

29 

28 

2 

12 

13 

1 

11 

30 

29 

4 

15 

16 

1 

12 

31 

30 

4 

16 

17 

1 

13 

32 

31 

6 

19 

20 

1 

14 

33 

32 

6 

20 

21 

1 

15 

34 

33 

8 

23 

24 

1 

16 

35 

34 

8 

24 

25 

1 

17 

27 

28 

54 

38 

37 

1 

18 

26 

29 

54 

39 

38 

1 

19 

29 

30 

56 

42 

41 

1 

20 

30 

31 

56 

43 

42 

1 

21 

31 

32 

58 

46 

45 

1 

22 

32 

33 

58 

47 

46 

1 

23 

33' 

34 

60 

50 

49 

1 

24 

34 

35 

60 

51 

50 

1 

25 

54 

53 

27 

36 

37 

1 

26 

55 

54 

29 

39 

40 

1 

27 

56 

55 

29 

40 

41 

1 

28 

57 

56 

31 

43 

44 

1 

29 

58 

57 

31 

44 

45 

1 

30 

59 

58 

33 

47 

48 

1 

F-gure  A. 4  Cantilever  Beam  Problem:  Input  for  the  5NT  Element 
Discretization. 


1  60 
2  61 
111 
1011 
2711 
3611 
5311 

1 

35200 
1250 
4150 
5000  5000 
1  9 

0.0 
53 


59  33  48  49 

60  35  51  52 


1 

3.1E-3 


100.0 


100.0 


Figure  A. 4 


Cantilever  Beam  Problem:  Input  for  the  5NT  Element 
Discretization  (Continued). 


STATIC  ANAL.  OF  A  CANT.  BEAM  USING  32  6NT  ELEMS .  6/20/85  BY  S.  MILLER 


85 

32 

1  5 

4000 

3.4E-5 

1  0 

2 

1 

32.2 

FT 

1 

3 

4  6 

-NODED  TRIANGULAR  CONTINUUM  ELEMENTS 

2.0 

2.0 

16.0 

4176 . 0E6 

0.0 

0.0 

0.0 

1 

0.0 

0.0 

17 

16.0 

0.0 

18 

0.0 

1.0 

34 

16.0 

1.0 

35 

0.0 

2.0 

51 

16.0 

2.0 

52 

0.0 

3.0 

68 

16.0 

3.0 

69 

0.0 

4.0 

85 

16.0 

4.0 

1 

1 

3  35 

2 

19 

18 

1 

4 

4 

13 

15  47 

14 

31 

30 

1 

4 

5 

3 

5  39 

4 

22 

21 

1 

4 

8 

15 

17  51 

16 

34 

33 

1 

4 

9 

3 

37  35 

20 

36 

19 

1 

4 

12 

15 

49  47 

32 

48 

31 

1 

4 

13 

3 

39  37 

21 

38 

20 

1 

4 

16 

15 

51  49 

33 

50 

32 

1 

4 

17 

35 

37  71 

36 

54 

53 

1 

4 

20 

47 

49  83 

48 

66 

65 

1 

4 

21 

37 

39  71 

38 

55 

54 

1 

4 

24 

49 

51  83 

50 

67 

66 

1- 

4 

25 

35 

71  69 

53 

70 

52 

1 

4 

28 

47 

83  81 

65 

82 

64 

1 

4 

29 

39 

73  71 

56 

72 

55 

1 

4 

32 

51 

85  83 

68 

84 

67 

1 

4 

111 

17 

-i 

6911 

17 

1 

1 

2 

51200 

1250 

6150. 

5000 

5000 

1 

17 

-1 

1 

3 

3. 

IE-3 

100.0 

1.0 

100 

69 


Figure  A. 5  Cantilever  Beam  Problem:  Input  for  the  6NT  Element 
Discretization. 
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APPENDIX  B 

Gauss-Legendre  Abscissae  and  Weight  Coefficients 


Table  B.l  shows  the  Gauss-Legendre  abscissae  and  weight 
coefficients  (±a  and  h)  which  are  used  in  the  Gaussian  quadrature  method 
of  numerical  integration.  These  coefficients  for  Gaussian  quadrature 
are  used  in  the  SAMS0N2  code  in  order  to  determine  the  lumped  mass 
matrix  [M],  the  internal  nodal  force  vector  < F1nt > .  the  internal  strain 
energy  of  the  system  U,  and  the  external  work  of  the  system  W  for  the 
8NQ  higher-order  isoparametric  continuum  element.  References  12,  14, 
and  15  contain  additional  information  about  Gaussian  quadrature  and  they 
also  contain  a  discussion  on  the  practical  applications  of  Gaussian 
quadrature  for  the  Isoparametric  element  formulation. 
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Table  B.l  Gauss-Legendre  Abscissae  and  Weight  Coefficients  for  a 
Particular  Order  of  Integration  n(n=l,10). 
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Reference  15. 
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APPENDIX  C 

nt  of  the  Elastic-Plastic  Stress-Strain  Matrix 


for  the  Plane  Strain  or  Axisymmetric  Material  Law 


The  general  equation  that  is  used  in  order  to  determine  the 

EP 

elastic-plastic  stress-strain  matrix  (C  )  for  biaxial  analysis  is  shown 
in  Equation  C.l  (see  pp.  388-389  of  Reference  14). 


where , 
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This  equation  was  formulated  under  the  assumptions  of  isothermal 
conditions.  Isotropic  hardening  and  an  associated  flow  rule.  The  values 
for  (q)  and  (p>  are  determined  using  the  expressions  shown  in 
Equation  C.2. 


(C.2) 


where , 


yield  function  that  specifies  the  state  of  multlaxial 
stress  which  corresponds  to  the  start  of  plastic  flow. 

normal  and  shearing  stress  components,  and 


plastic  strain  Increments. 


The  Ctr  matrix  for  the  biaxial  elastic-plastic  plane  strain  or 
axisynmetric  material  law  was  developed  by  substituting  the  expressions 
shown  in  Equation  C.3  into  Equations  C.l  and  C.2. 
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where  *  1  for  the  normal  stresses  and 
1J  *  0  for  the  shear  stresses,  and 
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The  development  of  CtH  was  performed  In  three  tasks.  The  first  task  was 
the  multiplication  of  CE  (Equation  4.39  in  Section  4. 2. 4. 2)  by  (q>.  The 
results  of  this  multiplication  are  shown  in  Equation  C.4. 


S11  4  T*  <s22  *  S33J 
$99  *  TT\T  +  S-«) 


CE(q> 


E(1  -  v) 


22  T  '*11  T  33' 

l-2v  - 
T^T  512 


S33  +  {S11  +  S22} 


(C.4) 


The  next  task  was  the  multiplication  of  (plT  by  (q)  with  the  result 


shown  in  Equation  C.5. 
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The  final  task  involved  the  substitution  of  the  results  shown  in 
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Equations  C.4  and  C.5  into  Equation  C.l.  This  task  resulted  in  the  C 
matrix  shown  in  Table  C.l.  This  matrix  was  used  in  subroutine  STRES3 
according  to  the  discussion  in  Sections  4. 2. 4. 2  and  4. 2. 4. 3. 
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